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I.  INTRODUCTION 

During  the  last  two  decades  the  growth  of  finite  element  software 
for  structure  applications  has  experienced  a  vast  changing  environment. 
In  the  1960's  research  efforts  in  finite  element  analysis  were  primarily 
concentrated  on  the  development  of  discretization  techniques  for  various 
structural  geometries.  In  this  connection,  numerious  special  purpose 
research  programs  were  written  for  testing  of  the  related  numerical 
algorithms.  Later,  efforts  were  evolved  into  the  development  of  large 
scale  general  purpose  computer  programs  such  as  NASTRAN,  STRUDL,  etc. 

In  the  1970's  attentions  were  turned  to  the  applications  of  the  finite 
element  method  to  nonlinear  problems.  These  include  understanding 
the  nature  of  various  nonlinear  systems,  investigations  of  nonlinear 
solution  algorithms  and  development  of  both  special  and  general  purpose 
programs  such  as  HONDO,  PLANS,  NONSAP  and  ADINA.  Since  nonlinearities 
in  a  structural  system  may  have  a  wide  range  of  meaning,  it  is  difficult 
to  write  a  software  in  complete  generality,  which  can  cover  all  the 
possible  nonlinear  features  that  a  structure  may  encounter.  Instead, 
most  finite  element  programs  were  written  with  a  limited  scope;  new 
features  and  solution  methods  are  added  as  they  become  available. 

Since  tremendous  amount  of  manpower  and  financial  resources  have 
already  been  invested  in  the  development  of  various  large  scale  general 
purpose  finite  element  software,  a  significant  cost  is  further  necessary 
in  order  to  maintain  and  upgrade  these  packages  coping  with  the  new 
computational  technology  and  changing  hardware  configurations.  Due 
to  the  high  cost  associated  with  the  program  maintenance  and  potential 
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budget  limitations,  there  is  an  apparent  need  to  conduct  a  coordinated 
effort  for  further  software  development  and  maintenance.  Consequently, 
several  user's  groups,  e.g.,  SAP,  STRUDL,  ADINA,  etc.,  have  been  formed 
for  information  and  cost  sharing.  Also  recognizing  such  a  need,  an 
ISEG  (Jnteragency  Software  Evaluation  Group)  associated  with  the  armed 
forces  and  various  government  agencies  was  formed  for  software  user's 
coordination.  The  immediate  goai  :f  ISEG  is  to  identify  application 
software  packages  and  perform^ in-depth  evaluation  of  these  packages. 
Initially,  such  effort  is  concentrated  on  the  structural  mechanics 
software  due  to  the  advanced  status  in  this  area.  The  selection  criteria 
and  evaluation  procedure  have  been  proposed  by  Nickel  1  [1].  As 

a  part  of  the  ISEG  effort,  a  nonlinear  finite  element  program  called 
ADINA  [2]  was  selected  for  the  software  evaluation. 

Among  many  nonlinear  finite  element  programs  that  are  available, 
ADINA  represents  one  of  the  several  which  have  been  used  extensively 
by  government  agencies  and  private  companies.  This  program  is  an  ex¬ 
tended  version  of  another  nonlinear  code,  called  NONSAP,  which  was 
developed  at  the  University  of  California,  Berkeley.  Later,  NONSAP 
was  extended  and  modified  by  Bathe  of  MIT,  one  of  the  original 
NONSAP  developers,  to  include  additional  features  such  as  out-of-core 
solution  scheme,  double  precision  arithmetic,  expand  element  and  ma¬ 
terial  libraries,  etc.  The  extended  version  was  then  renamed  as  ADINA. 
(Automatic  Dynamic  Incremental  Nonlinear  Analysis). 

Use  of  nonlinear  finite  element  computer  programs  can  be  a  difficult 
task  due  to  the  numerical  complexity  of  nonlinear  analysis  in  general 
and  the  limitation  of  the  program  in  particular.  When  an  analysis 


T 


is  performed,  computation  may  be  interrupted  for  a  number  of  reasons 
and  it  is  hard  to  pin-point  the  exact  cause.  In  many  cases,  unsuccess¬ 
ful  runs  were  due  to  the  mis-use  of  the  program  by  the  analyst  resulting 
from  inadequate  experience  or  the  lack  of  proper  understanding  the 
limitations  of  the  code.  Any  supportive  information  leading  to  the 
effective  usage  of  the  nonlinear  program  would  be  helpful  from  the 
user's  standpoint.  It  is  therefore  the  intent  of  this  research  effort 
to  conduct  a  systematic  review  and  evaluation  of  ADINA  in  terms  of 
its  capability  and  limitations  for  solving  nonlinear  structural  problems. 

By  following  the  procedure  outlined  in  [1],  the  evaluation  work 
consists  of  the  following  items: 

1)  Theoretical  basis,  assumptions  and  numerical  approaches 
adopted  in  the  program. 

2)  Program  architecture  in  terms  of  program  design,  use  of 
control  variables  and  data  structure. 

3)  Review  of  nonlinear  static  and  dynamic  solution  algorithms. 

4)  Identifications  of  typical  algorithmic  limitations  encountered 
for  both  the  nonlinear  static  and  dynamic  analyses. 

Comments  for  the  first  two  items  are  enumerated  in  Part  I  of  this  report, 
whereas  the  results  obtained  for  the  last  two  items  are  discussed  in 
Part  II. 

Since  at  the  beginning  of  this  work,  the  1977  ADINA  was  made  availa¬ 
ble  by  the  developer  to  the  investigators.  Therefore,  all  the  comments 
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and  results  are  limited  to  this  version  of  ADINA.  Although  a  more 
recent  version,  1978  ADINA*,  was  released,  it  was  not  made  available 
for  carrying  out  this  study.  Nevertheless,  a  description  on  the  1978 
ADINA  is  outlined  in  section  II  of  this  report  based  on  the  revised 
user's  manual . 


*  The  1978  ADINA  was  released  to  the  ADINA  user  group  in  the  later 


part  of  1979. 


II.  GENERAL  DESCRIPTION 


ADINA  [2,3]  is  a  general  purpose  finite  element  program  which 
was  developed  for  conducting  both  linear  and  nonlinear,  static  and 
dynamic  analysis  for  a  range  of  structures.  The  program  was  extended 
from  another  nonlinear  structural  analysis  program,  i.e.  its  predecessor 
NONSAP  [4],  which  was  released  to  the  public  domain  in  1974. 

Although  ADINA  can  be  used  for  both  linear  and  nonlinear  struc¬ 
tural  analysis,  the  major  function  of  the  program  is  its  unique  capa¬ 
bility  in  performing  nonlinear  analysis.  Therefore,  it  is  viewed 
primarily  as  a  nonlinear  code  compared  to  other  linear  codes  such  as 
NASTRAN  [5],  STRUDL  [6]  and  SAP  [7]. 

The  analysis  capability  of  a  general  purpose  finite  element  program 
may  be  best  described  by  examining  four  different  areas,  namely 

1)  Analysis  procedures 

2)  Element  library 

3)  Material  library 

4)  Other  features 

The  above  items  for  the  1977  and  1978  versions  of  ADINA  in  comparison 
with  NONSAP  are  summarized  in  Table  1.  Each  of  these  areas  is  dis¬ 
cussed  in  the  following  sections  in  accordance  with  the  ADINA  user's 
manual  [2]. 
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1 .  Analysis  Procedures 

In  addition  to  its  linear  analysis  capability,  the  analysis 
procedures  consist  of  three  major  areas: 

1)  Nonlinear  static  analysis 

2)  Nonlinear  dynamic  analysis 

3)  Frequency  analysis 

Nonlinearities  include  both  nonlinear  materials  and  large  deforma¬ 
tions.  However,  the  large  deformation  that  can  be  utilized  is 
limited  to  large  displacement  (or  rotation)  but  small  strain. 

This  is  due  to  the  fact  that  the  stress  and  strain  measures  adopted 
in  the  constitutive  formulation  of  the  code  are  not  in  compliance 
with  the  requirements  of  large  strain  analysis.  A  further  discus¬ 
sion  of  this  point  will  be  given  in  section  III. 2.  The  solution 
method  for  solving  nonlinear  problems  (both  static  and  dynamic) 
is  based  on  the  incremental  approach  together  with  equilibrium 
iterations.  For  each  load  (or  time)  increment,  tangent  stiffness 
is  used  for  solving  the  displacement  equilibrium  equations  and 
the  stiffness  of  the  structure  remains  unchanged  during  iterations. 
Therefore,  the  method  is  equivalent  to  the  well-known  modified 
Newton-Raphson  scheme  for  solving  nonlinear  equations  by  numerical 
approach  [8]. 

For  dynamic  analysis,  the  equilibrium  equations  of  the  finite 
element  system  are  solved  by  step-by-step  numerical  integration 
in  time,  for  which  both  the  implicit  and  explicit  integration 
(or  central  difference)  methods  are  adopted  in  the  program. 


V.  it'jf  AkL 
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Table  1  Summary  of  Analysis  Options  and 
Comparison  with  NONSAP 


1978 

1977 

ADINA 

ADINA 

NONSAP 

Analysis  Procedures 

Nonlinear  Static  Analysis 

X 

X 

X 

Nonlinear  Dynamic  Analysis: 

Implicit  integration  method 

X 

X 

X 

Explicit  integration  method 

X 

X 

Frequency  Analysis 

X 

X 

X 

Element  and  Material  Library 

Truss  Element: 

X 

X 

X 

Linear  elastic 

X 

X 

X 

Nonlinear  elastic 

X 

X 

X 

Thermal -elastic 

X 

X 

Elastic-plastic 

X 

X 

Thermal -elastic -pi  as tic-creep 

X 

X 

2/D  Continuum 

X 

X 

X 

Isotropic  linear  elastic 

X 

X 

X 

Orthotropic  linear  elastic 

X 

X 

X 

Isotropic  thermal -el astic 

X 

X 

Curve  description  model 

X 

X 

X 

Concrete  model 

X 

X 

Elastic-plastic  Drucker-Prager 

X 

X 

X 

Elastic-plastic  Von  Mises 

X 

X 

X 

Thermal -el  as tic-pl as tic -creep 

X 

X 

Mooney-Rivl in  material 

X 

X 

X 

3/D  Continuum 

X 

X 

X 

Isotropic  linear  elastic 

X 

X 

X 

Orthotropic  linear  elastic 

X 

X 

Isotropic  thermal-elastic 

X 

X 

Curve  description  model 

X 

X 

X 

Concrete  model 

X 

X 

Elastic-plastic  Drucker-Prager 
Elastic-plastic  Von  Mises 

X 

X 

X 

Thermal -elastic-plastic -creep 

X 

X 
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1978 

ADINA 

1977 

ADINA 

NONSAP 

Element  and  Material  Library 

Beam  Element 

X 

X 

Linear  elastic 

X 

X 

Elastic-plastic 

X 

X 

Shell  Element 

Linear  elastic 

Elastic-plastic 

2/D  Fluid 

X 

X 

X 

X 

(Invicid  compressible  fluid 
with  constant  bulk  modulus) 

3/0  Fluid 

X 

(Invicid  Compressible  fluid 
with  constant  bulk  modulus) 

Other  Features 

Out-of-core  linear  equation  solver 

X 

X 

Analysis  restart 

X 

X 

X 

Data  check 

X 

X 

X 

Nodal  force  Loading 

X 

X 

X 

Pressure  loading* 

X 

X 

Thermal  loading 

X 

X 

Displacement  boundary  Conditions 

X 

X 

Companion  thermal  analysis  code 

X 

X 

Data  Card  image  printout 

X 

X 

Substructuring  X 


*8^6^  of  geometric  change  was  not  included. 


In  the  implicit  integration,  either  Nemark's  g-method  or  Wilson 
0-method  can  be  applied.  Options  are  available  to  generate  either 
the  lumped  mass  or  consistent  mass  matrix  for  the  structure  in 
the  same  way  as  a  linear  problem.  However,  if  the  central  differ¬ 
ence  method  is  used,  only  the  lumped  mass  matrix  is  permitted. 

In  addition,  concentrated  masses  can  be  added  at  selected  degrees 
of  freedom.  To  consider  the  damping  effect,  only  concentrated 
nodal  dampers  are  included  and  this  information  must  be  specified 
for  each  direction  of  those  nodes  having  the  damping  effect. 

In  the  frequency  analysis,  two  solution  algorithms,  namely 
the  determinant  search  method  and  the  subspace  iteration  technique, 
are  available.  From  the  analysis,  the  lowest  p  eigenvalue  and 
corresponding  mode  shapes  are  calculated.  An  over-relaxation 
scheme  with  a  shifting  strategy  is  also  available  in  the  subspace 
iteration  technique  to  accelerate  the  convergence  of  the  solution. 
Since  the  eigen-solution  in  ADINA  is  performed  outside  the  time- 
integration  loop,  the  stiffness  matrix  used  in  the  analysis  is' 
essentially  linearly  elastic  even  though  the  problem  in  question 
may  be  nonlinear.  The  mass  matrix  can  be  either  lumped  or  consis¬ 
tent.  No  damping  effect  is  considered  in  the  frequency  analysis. 

As  pointed  out  in  the  ADINA  manual,  the  frequency  analysis 
is  desirable  for  estimating  the  resonance  conditions  of  a  vibrating 
system,  or  for  selecting  a  suitable  time  step  At  in  a  direct  step- 
by-step  numerical  integration.  However,  such  estimation  is  valid 
only  for  a  linear  system. 
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2.  Element  Library 

Element  type  is  used  to  represent  the  geometric  configuration 
of  a  specific  structure.  Therefore,  the  structural  geometry  that 
can  be  handled  by  a  program  is  determined  by  the  types  of  element 
available  in  the  code.  In  the  1977  version  of  ADINA,  there 
are  basically  four  types  of  elements: 

1)  Three-dimensional  truss 

2)  Three-dimensional  beam 

3)  Two-dimensional  (2-D)  isoparametric  solid  (for  plane 
stress,  plane  strain  and  axisymmetric  deformations) 
with  4  to  8  nodes. 

4)  Three-dimensional  (3-D)  isoparametric  solid  with  8  to 
20  nodes. 

For  the  truss  and  beam  elements,  large  displacement  with  small 
strain  was  assumed  and  the  element  stiffnesses  were  derived  on  the 
basis  of  the  updated  Lagrangian  formulation  (UL)*.  Whereas  for  the 
2-D  and  3-D  solid  elements,  both  the  total  Lagrangian  (TL)  and' 
updated  Lagrangian  formulations*  are  available.  Although  a  structure 
can  be  modeled  by  a  combination  of  different  element  types,  cautions 
must  be  given  in  two  accounts: 

1)  In  the  case  of  large  deformation  analysis,  either  the 
total  or  the  updated  Lagrangian  formulation  should  be 
used  throughout  for  all  element  types. 

2)  If  a  structure  is  to  be  modeled  by  a  mixture  of  different 
element  types,  it  appears  that  the  user  can  only  combine 


*  Definitions  of  total  and  updated  Lagrangian  formulations  are  given 
in  Section  III. 
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the  truss  and  2-D  solid  elements.  Any  other  combination 
of  different  element  types  will  not  satisfy  the  compatibil¬ 
ity  condition  at  the  junction  of  the  nodes  due  to  the  dif¬ 
ference  in  degrees  of  freedom.  In  order  to  remove  this 
difficulty,  appropriate  constrain  conditions  must 
be  imposed  at  the  nodes  involved. 

3.  Material  Models 

In  finite  element  analysis,  the  definition  of  constitutive 
relations  of  a  material  is  needed  for  calculation  of  the  stiffness 
matrices  and  stress  components  at  the  element  level.  The  material 
models  included  in  ADINA,  as  seen  in  Table  1,  are  element  -  depen¬ 
dent.  That  is  to  say,  if  the  material  model  is  available  to  one 
element  type,  it  may  not  be  available  to  other  element  types. 

The  2/D  continuum  element  has  the  most  complete  list  of  material 
models  ranging  from  linear  to  nonlinear  stress-strain  relation¬ 
ships.  For  a  quick  reference,  the  material  models  as  opposed 
to  the  element  library  in  ADINA  are  shown  in  Table  2.  As  seen 
from  this  table,  the  element  that  has  the  most  material  model 
coverage  is  the  2/D  continuum. 

In  addition,  an  "element  birth  and  death"  option  is  available 
for  all  elements  and  material  models  in  ADINA.  The  element  birth 
option  will  activate  an  element,  which  is  initially  not  active, 
at  a  specified  time  of  birth.  On  the  other  hand,  the  death  option 
makes  an  active  element  to  become  inactive  at  its  time  of  death. 


Table  2.  Availability  of  Material  Models  vs. 


Element  Library 


Material  Models  Truss 

Isotropic  linear  elastic  x 

Orthotropic  linear  elastic 
Isotropic  thermo-elastic  x 

Nonlinear  elastic  x 

Curve  description  model 
Concrete  model 


Drucker-Prager  model + 

von  Mises  elastic-plastic  with  x 
isotropic  or  kinematic 
hardening 

Thermo-elastic-plastic  and  x 

creep  with  isotropic  or 
kinematic  hardening 

Mooney-Rivlin  material 


2/D  Continuum 


3/D  Continuum 
* 


Beam 

x 


x 


x 


x 


x 


x 

x 

x 

x 


x 

x 

x  x 


x  x 

x+ 


Notes: 

*  In  the  case  of  large  deformation  analysis,  the  isotropic  linear 
elastic  model  can  be  used  for  both  TL  and  UL  formulations,  whereas 
other  material  models  are  limited  to  TL  formulation  only. 

+  Drucker-Prager  model  can  only  be  used  for  two-dimensional  plane 
strain  or  axisymmetric  deformation  analysis 

t  Mooney-Rivlin  material  can  only  be  used  for  plane  stress  analysis 
with  the  TL  formulation. 
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According  to  the  manual,  these  options  are  useful  in  the  analysis 
of  construction  or  excavation  processes.  It  is  not  clear,  however, 
to  the  evaluators  what  theoretical  implication  these  options  will 
cause.  Furthermore,  sudden  declaration  of  element  birth  or  death 
during  the  critical  stage  of  loading  may  introduce  numerical  in  - 
stability  when  iterations  are  being  performed.  Therefore,  the 
user  must  apply  the  element  birth  or  death  option  with  some  care. 

4.  Other  Features 


As  seen  in  Table  1,  one  of  the  major  improvements  of  ADINA 
over  its  predecessor  NONSAP  is  the  availability  of  the  out-of-core 
linear  equation  solver.  This  option  provides  greater  flexibility 
for  the  user  to  handle  large  size  application  problems  on  most  of 
the  main  frame  computers.  Other  useful  improvements  over  the 

NON SAP  include 

1)  Calculation  of  pressure  loading 

2)  Thermal  loading 

3)  A  companion  temperature  analysis  code  ADINAT  [9] 

From  the  user's  standpoint,  the  biggest  drawback  of  ADINA 

is  probably  the  lack  of  a  comprehensive  pre-  and  post-processing 
capability.  The  pre-processing  in  terms  of  automatic  mesh  genera¬ 
tion  and  on-line  display  of  finite  element  grid  with  editing  capa¬ 
bility  is  extremely  useful  for  conducting  nonlinear  analysis 
of  large  scale  structures.  Nevertheless,  the  program  provides 
the  option  to  read  the  input  information  from  a  porthole  tape 
on  which  the  data  can  be  generated  from  an  independent  preprocessor 


Data  check  option  is  available,  but  it  is  far  from  sufficient. 
The  current  option  only  prints  out  the  input  data.  Some  prelimi¬ 
nary  calculations,  e.g.  the  area  or  volume  of  the  element,  deter¬ 
minant  of  the  Jacobian  matrix,  etc.,  should  be  included  to  check 
any  mistakes  made  in  nodal  coordinates  and  element-nodal  number 
relationship.  To  further  demonstrate  this  point,  for  example, 
the  most  frequently  encountered  error  message  in  using  ADINA  is 
the  one  printed  in  the  equation  solver  "COLSOL",  i.e. 

"STIFFNESS  NOT  POSITIVE  DEFINITE 
NONPOSITIVE  PIVOT  FOR  EQUATION  XX 
PIVOT  =  (Value  of  diagonal  term)" 

The  reason  for  causing  this  error  can  be  many.  Typically,  they 
are 

1)  Mistakes  in  nodal  coordinates 

2)  Mistakes  in  element  -  nodal  number  definitions 

3)  Erroneous  material  constants 

4)  Accumulated  numerical  errors  due  to  insufficient  iterations 

5)  The  problem  has  reached  structural  or  material  instability. 
Because  of  the  above  mentioned  sources  of  errors  which  are  diffi¬ 
cult  to  pin-point,  it  is  imperative  for  a  nonlinear  analysis  to 
go  through  a  more  thorough  data  checking  process. 

5.  Review  of  Documentations 

ADINA  documents  consist  of  primarily  the  program's  instruc¬ 
tion  manual  [2]  and  the  theoretical  manual  [3].  In  the  theoretical 
manual,  a  description  of  the  theoretical  basis,  finite  element 


15 


discretization,  definitions  of  material  models,  numerical  approach 
and  solution  of  sample  problems  are  systematically  presented. 

This  manual  provides  a  sound  basis  for  the  numerical  algorithm 
undertaken  in  the  code.  The  instructional  manual  consists  of 
three  main  parts:  i)  Description  of  ADINA,  ii)  Sample  problems, 
and  iii)  User's  instruction.  These  documents  have  been  reviewed 
and  several  conments  are  drawn. 

Commentary: 

1)  Theoretical  manual  is  well  documented,  in  which  the 
element  stiffness  matrices  and  definitions  of  various  material 
models  adopted  in  the  code  are  available  for  users  who  wish  to 
understand  the  coding  procedure. 

2)  For  large  deformation  analysis,  both  .the  total  Lagrangian 
(TL)  and  the  updated  Lagrangian  (UL)  formulations  are  given  for 
2/D  as  well  as  3/D  continuum  elements.  At  several  places  in  the 
manual,  it  was  stated  that  the  choice  between  the  TL  and  UL  formu¬ 
lations  lies  essentially  on  their  relative  numerical  efficiency. 
From  the  stiffness  formulation  of  the  2/D  element,  for  example, 

the  linear  part  of  the  strain-displacement  transformation  matrix 
for  the  TL  formulation  is  a  full  matrix,  whereas  for  the  UL  formu¬ 
lation  is  sparse.  Consequently,  the  UL  formulation  appears  to 
be  more  numerically  efficient  than  the  TL  formulation.  Following 
this  implication,  the  user  may  become  pursuaded  to  use  the  UL 
formulation  for  all  situations  without  considering  other  impor¬ 
tant  factors  such  as  the  suitability  of  constitutive  relations 
for  a  given  material,  etc. 


3)  The  user's  instruction  manual  was  clearly  written  and 
fairly  easy  to  follow.  Even  for  a  new  user  of  the  ADINA  code, 
the  input  data  can  be  prepared  according  to  the  instructions  with¬ 
out  much  confusion. 

4)  For  the  data  input  some  of  the  control  variables  used 
in  the  program  could  be  consolidated  so  that  the  number  of  input 
blank  cards  can  be  reduced.  For  example  if  a  control  variable 

is  provided  on  input  Card  1  of  Section  II  to  distinguish  the  static, 
dynamic  or  frequency  analysis  option,  then  the  blank  cards  #2, 

#3,  and  #4  can  be  eliminated  when  a  static  analysis  is  requested. 

5)  Although  ADINA  does  not  have  a  sophisticated  pre-  or  post¬ 
processor  package,  the  program  has  already  set  up  the  ground  work 

to  read  the  nodal  and  element  data  from  a  tape  created  by  a  pre¬ 
processor.  Also,  the  analysis  results,  such  as  nodal  displace¬ 
ments  and  element  stresses,  generated  by  ADINA  can  be  written 
on  an  output  tape  for  plotting  or  further  evaluation  by  a  post¬ 
processor.  Therefore,  the  program  can  easily  be  interfaced  with 
a  pre-  and  post-processing  software. 

6)  In  an  attempt  to  achieve  better  computational  efficiency, 
the  program  provides  the  options  to  the  users  that  stiffness  re¬ 
formation  and  equilibrium  iterations  can  be  requested  for  certain 
blocks  of  time  steps.  However,  it  is  not  clear  how  these  options 
can  be  effectively  utilized  unless  the  analyst  knows  the  conver¬ 
gence  characteristics  of  a  nonlinear  problem  a  priori.  The  re¬ 
viewers  feel  that  the  option  on  iteration  time  blocks  may  have 
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greater  application  in  handling  the  elastic-plastic  problems  for 
which  convergent  solution  is  difficult  to  obtain  when  the  struc¬ 
ture  is  experiencing  plastic  loading  and  then  unloading.  This 
point  will  be  further  demonstrated  in  Part  II  of  this  report  [10]. 

III.  THEORETICAL  BASIS 

As  mentioned  earlier,  the  ADINA  program  can  perform  both  linear 
and  nonlinear,  static  and  dynamic  analysis  for  a  range  of  structures. 
Both  geometric  and  material  nonlinearities  were  included  in  its  theo¬ 
retical  formulations.  In  this  section,  the  basic  theory,  constitutive 
relations,  finite  element  idealization,  and  numerical  scheme  for  both 
static  and  dynamic  analysis  will  be  reviewed  and  commented. 

1 .  Theory 

The  solution  approach  adopted  in  ADINA  is  basically  an  in¬ 
cremental  method  for  solving  nonlinear  problems,  whether  static 
or  dynamic.  For  each  loading  increment,  a  tangent  stiffness  is 
used  in  the  solution  of  displacement  finite  element  equations 
and  the  solution  is  improved  by  invoking  equilibrium  iterations, 
while  the  stiffness  is  kept  constant.  This  method  corresponds 
to  the  so-called  modified  Newton-Raphson  method.  For  the  treat¬ 
ment  of  large  deformation,  both  the  total  Lagrangian  and  updated 
Lagrangian  formulations  are  adopted.  A  brief  review  on  each  of 
the  two  formulations  is  given  as  below. 
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The  TL  formulation  for  large  deformation  finite  element  analysis 
has  been  adopted  by  a  number  of  researchers,  and  the  works  that  can 
be  cited  are,  for  example,  Martin  [11],  Oden  [12],  Hibbitt,  Marcal 
and  Rice  [13],  McNamara  [14],  Sharif i  and  Yates  [15],  Needleman  [16] 
and  Hutchinson  [17]  also  used  the  TL  formulation  with  convected  coordi¬ 
nates  to  study  finite  strain  plasticity  problems.  For  the  purpose 
of  discussion,  a  brief  review  of  this  formulation  is  given. 

For  a  nonlinear  finite  element  system,  the  incremental  equations 
of  equilibrium  are  generally  derived  from  a  virtual  work  principle. 

In  reference  to  a  total  Lagrangian  description,  the  equilibrium  con¬ 
dition  of  a  structure  at  time  t+At  can  be  stated  by  the  principle  of 
virtual  displacement  [2,18],  i.e. 


f  C..  e  6  e .  .  °d v  +  f  S--  6  n ■  •  °d i 

J  o  ijrs  o  rs  o  ij  J  o  ij  o  'ij 

0\/  0\/ 


_  t+AtD  f  tC  r  Q  0 

R  '  i  oSij  5oeij 


d  v 


(1) 


where  t+AtR 


=  1  6UK  °da  +  J  °p  t  +  AofK  6  u K  °dv  (2) 
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C.  . 

o  ijrs 


t+At. 

0  1 


=  Area  of  body  in  the  configuration  at  time  o. 

=  Component  of  tangent  constitutive  tensor  at  time  t  referred 
to  the  configuration  at  time  o. 

=  Incremental  Green-Lagrange  strain  tensor  referred  to  the 
configuration  at  time  o 

j 

=  2"  P iola-Kirchhoff  stress  tensor  in  the  configuration 

at  time  t 

=  Linear  part  of  strain  increment  e.. 

o  1J 

=  Nonlinear  part  of  strain  increment  e. • 

o  i  j 

=  Component  of  surface  traction  vector  per  unit  area,  at 
time  t+At,  referred  to  the  configuration  at  time  o 


ui  =  Component  of  incremental  displacement 


t+Atf 
o  K 


p  =  Specific  mass  of  the  body  in  the  configuration  at  time  o 

^  =  Component  of  body  force  vector  per  unit  mass,  at  time 

t+At,  referred  to  the  configuration  at  time  o. 


In  the  case  of  dynamic  analysis,  a  term  due  to  inertia,  i.e.  {-/°pt+AtuK 
6uk  °dv),  must  be  added  to  the  left  hand  side  of  Eq.  (1).  The  stress  and 
strain  components  at  time  t+At  are  calculated  from  the  following  relations: 


t+Ats. .  =  ts. .  +  s. . 

o  ij  o  1J  o  1J 
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t+Atc  =  V.  +  e.  .  «> 


obi  j 

E  .  .  + 

0  IJ 

oei  j 

\  ^  / 

and 

oei  j 

=  e. .  + 

O  IJ 

onij 

(5) 

oei  j 

■  Kyj 

ouj,i  +  o  uK,i  ouK,j  +  ouK,i  o  uK,j^ 

(6) 

onij 

1 

2  oUK,i 

oUK,j 

(7) 

where 

u .  - 
0  1  ,J 

_  Vi 

Vj 

(8) 

t 

oui ,  j 

3*u. 

O  1 

3  x. 
o  J 

(9) 

In  order  to  make  Eq.  (1)  solvable,  the  increment  of  strain  components 


is  approximated  by 


e .  . 
o  ij 


"  oeij 


(10) 


Therefore,  Eq.  (1)  becomes 


[  C..  e  e..  °dv  +  [  V.  <5  n--  °dv 
}Q  ijrs  o  rs  o  ij  J  o  ij  o  ij 

_  t+AtD  [  tS. .  6  e . .  °dv 

=  R  -  i  o  ij  o  ij 


(11) 


This  is  the  virtual  work  equation  based  on  which  the  total  Lagrangian 
finite  element  formulation  was  derived  in  ADINA. 

It  is  seen  that  the  major  assumption  made  in  arriving  at  Eq.  (11)  is 
the  linearization  of  the  incremented  strain  tensor  by  replacing  j  by 
Qe^j  in  the  first  term  on  the  left  hand  side  of  Eq.  (1).  The  implication 
of  this  approximation  is  that  the  strain  increment  in  an  element  of 
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material  for  each  time  (or  loading)  step  must  be  restricted  to  be  small. 
Otherwise,  the  equilibrium  condition  of  the  structure  cannot  be  satis¬ 
fied.  Another  approximation  made  is  the  evaluation  of  virtual  exter¬ 
nal  work  on  the  right  hand  side  of  Eq.  (1).  Generally,  in  the  case 
of  large  deformation  the  virtual  external  work  should  consist  of  two 
parts:  a  part  due  to  the  change  of  load  value  itself,  and  another 
part  due  to  the  change  in  geometry  while  the  load  value  remains  con¬ 
stant  [13,  19,  20].  In  the  ADINA  program,  the  latter  part,  i.e.,  the 
effect  of  deformation  dependent  load,  was  neglected. 

The  updated  Lagrangian  formulation  uses  the  current  conf iguration 
of  a  body,  i.e.  at  time  t,  under  consideration  as  a  reference  state 
to  obtain  the  state  of  deformation  at  next  time  step  t+At.  This  ap¬ 
proach  has  previously  been  used  by  Yaghmai  and  Popov  [21],  Sharif i 
and  Popov  [22],  Hofmeister,  Greenbaum  and  Evenson  [23].  Osias  and 
Swedlow  [24]  applied  the  UL  formulation  for  the  solution  of  large 
elastic-plastic  problems  with  a  Galerkin  procedure  which  admits  non- 
syrrmetric  constitutive  relations.  Based  on  Hill’s  work  [25], 

McMeeking  and  Rice  [26],  Yamada  [27],  Yamada,  Hirahawa  and  Wifi  [20], 
and  Namat-Nasser  and  Taya  [28]  have  also  presented  a  UL  formulation 
in  dealing  with  large  strain  plasticity  problems.  An  improved  UL 
formulation  to  deal  with  the  numerical  difficulty  arose  due  to  the 
incompressible  behavior  in  large  strain  plasticity  was  presented  by 
Atluri  [29],  Murakawa  and  Atluri  [30]. 

According  to  ADINA,  the  virtual  work  principle  corresponding  to 
the  UL  formulation  is  given  by 


jL.  *  .  _ M:>  »•  _ 
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J  tC. . 

ty  ^ 

ters  teij  ‘d*  +  I  Sj  *tni j  *dv 

Lv 

=  t+At 

R  •  1  Sj  steij  dv 
cv 

(12) 

where 

tCijrs 

=  Component  of  constitutive  tensor  at  time  t 

referred 

to  the  configuration  at  time  t. 

t£ij 

=  Incremental  Green-Lagrangian  strain  tensor 

referred 

to  the  congi figuration  at  time  t 

teij 

=  Linear  part  of 

L  1 J 

tnij 

=  Nonlinear  part  of  .  e . . 

L  1 J 

IJ 

=  Cauchy  stress  tensor  in  the  configuration  at  time  t. 

tV 

=  Volume  of  the  body  occupied  at  time  t. 

The  stress  and  strain  components  at  time  t+At  are  calculated  from 

t+At,- 
t  ij 

t  .  c 

ij  tSij 

(13) 

t+At 
t  ij 

te i j  "  tei j  +  tnij 

(14) 

tei  j 

2  ^tui,j  +  tuj  ,i ^ 

(15) 

tnij  = 

2  tUK,i  tUK,j 

(16) 

with 

tui,j  = 

aui 

(17) 

▲ 
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Therefore,  Eq.  (12)  forms  the  basis  for  the  UL  finite  element  formula¬ 
tion  in  ADINA.  Similar  to  the  TL  formulation,  this  equation  involves 
two  approximations:  i)  the  strain  increment  e..  in  the  first  term 

t  1  J 

on  the  left  hand  side  of  Eq.  (12)  was  linearized  to  become  .e..,  and 

t  1 J 

ii)  the  virtual  variation  of  external  work  due  to  any  geometric  change 
was  neglected. 

It  is  also  important  to  note  that  for  large  deformation  analysis 
appropriate  definition  for  the  constitutive  tensor  tCijrs  must  be 
assigned  to  relate  the  stress  .S. .  and  strain  .e..,  and  this  point 

t  1  J  t  I J 

will  be  further  discussed  in  the  next  section. 

Although  both  the  TL  and  UL  formulations  were  presented  in  the 
ADINA  theoretical  manual  [2],  they  are  not  available  for  all  the  ele¬ 
ments.  As  seen  in  Table  3,  only  the  2/D  and  3/D  continuum  elements 
were  derived  for  both  formulations. 
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Table  3.  Large  Deformation  Formulations 


Element  Type 
Truss  Element 
2/D  Continuum 
3/D  Continuum 
Beam  Element 
Shell  Element* 
2/D  Fluid* 

3/D  Fluid* 


Total  Lagrangian 

X 

X 

X 


Updated  Lagrangian 

X* 

X+ 

X  + 

X* 

X 

X 


Notes : 

*  .  Large  displacements  but  small  strains. 

+  .  Only  isotropic  linear  elastic  material  model  can  be  used, 
t  .  Element  is  available  in  1978  version. 


v. 


Mtfuki 
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2.  Constitutive  Relations 

One  of  the  important  features  of  a  finite  element  program  is  its 
ability  to  model  the  responses  of  various  materials.  In  this  connection, 
ADINA  has  included  a  number  of  material  models  ranging  from  linear  to 
nonlinear  behavior.  With  the  consideration  of  the  kinematics  in  its 
general  form  as  seen  in  the  virtual  work  equations,  i.e.  Eqs.  (11)  and 
(12),  it  appears  that  the  program  can  be  applied  to  any  type  of  nonlinear 
analysis  of  structural  problems  involving  either  nonlinear  material, 
large  deformation  or  both.  However,  care  must  be  given  when  both  material 
and  geometric  nonl inearities  are  to  be  considered  simultaneously.  In 
fact,  the  version  of  ADINA  being  reviewed  is  only  suitable  for  large 
deformation  analysis  in  the  sense  of  large  displacement  but  small  strain 
(except  for  the  Mooney-Rivl in  material  model).  The  discussion  on  this 
point  can  be  made  by  reviewing  some  of  the  fundamentals  in  constitutive 
theory  as  follows. 

The  constitutive  relations  of  a  material  can  be  expressed  in 
different  forms  in  conjunction  with  the  stress  and  strain  measures  chosen. 
Following  ADINA's  notations,  three  different  expressions  are  to  be 
discussed:  Incremental  Piola-Kirchhof  stress  vs.  incremental  Green  strain, 
updated  incremental  Piola-Kirchhof  stress  vs.  updated  incremental  Green 
strain,  and  Jaumann  rate  of  Kirchoff  stress  vs.  deformation  rate.  These 
relationships  are  given  by 


where 


Jaumann  (frame  invariant)  rate  of  Kirchhoff  stress 
tensor  in  the  configuration  at  time  t. 


In  the  above  equations,  the  consitutive  tensors  are  related  through 


C  =  _P  0 

o  i jrs  t  t  xi ,m 

P 

°x  °x  °x  C 

txj,n  t  r,p  tXs»q  t mnpq 

(22) 

and 

c  =  cv 

t  ijrs  t  ijrs 

T .  5  ■  -  -  5 •  +  .  .  6 

lr  ujs  jr  is  tj  rs 

(23) 

where  t  =  Specific  mass  of  the  body  in  the  configuration  at  time  t 

t°x.  =  Derivative  of  coordinate  in  the  configuration  at  time  o  with 
1,m  respect  to  the  coordinate  at  time  t 

“  T -  (  xi > 

3xm 

m 

*t..  =  Cauchy  stress  tensor  at  time  t. 

*  w 

The  choise  on  the  use  of  Eq.  (18),  (19)  or  (20)  depends  on  the  type  of 
material  in  question.  For  example  if  a  material  is  hypereleastic  for 
which  the  constitute  relations  are  derivable  from  a  strain  energy  function 


then  Eq.  (18)  would  be  the  natural  choice  for  the  description  of  material 
response.  On  the  other  hand,  for  hypoelastic  materials  such  as  metal 
plasticity,  the  constitutive  behavior  depends  on  the  stress  and  strain 
history  developed  in  the  structure  and  the  natural  description  of  the 
material  response  would  be  to  relate  the  frame  invariant  stress  rate 
and  the  deformation  rate  as  given  in  Eq.  (20).  As  discussed  by  Hibbitt, 
Marcal  and  Rice  [13]  in  accordance  with  Hill's  interpretation  [27], 

Eq.  (20)  reduces  to  the  true  stress  vs.  logarithmic  strain  relation 
in  simple  tension  of  metals. 

From  the  theoretical  standpoint,  it  is  indifferent,  in  the  finite 
element  derivation,  to  choose  one  form  of  constitutive  description 
over  the  other,  as  long  as  the  transformation  relations,  i.e.  Eq.  (22) 
or  (23)  are  properly  defined  in  the  program.  However,  from  the  numeri¬ 
cal  point  of  view,  it  is  preferable  to  use  the  stress  and  strain  mea¬ 
sures  in  the  finite  element  computations  same  as  those  for  the  consti¬ 
tutive  description.  Although  the  transformations  between  the  consti¬ 
tutive  tensor  can  be  made  according  to  Eq.  (22)  or  (23),  such  calcula¬ 
tions  are  very  time  consuming  and  this  point  was  previously  discussed 
in  references  [13,  27].  Therefore,  for  large  strain  analysis  of  elas¬ 
tic  materials,  it  is  more  efficient  to  use  the  TL  approach;  whereas 
for  large  strain  elastic-plastic  analysis,  the  l)L  approach  would  be 
the  natural  choice. 

As  we  turn  our  attention  to  the  ADINA  program,  most  of  the  nonlinear 
material  models  are  restricted  to  the  TL  formulation  and  type  of  trans¬ 
formation  between  the  constitutive  tensors  as  given  in  Eqs.  (22)  and  (23) 
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is  not  available  in  the  code.  It  is  therefore  concluded  that  the  present 
version  of  ADINA  is  only  suitable  for  small  strain,  nonlinear  material 
analysis  or,  at  most,  for  large  displacement  nonlinear  material  analysis. 

The  second  point  worth  noting  is  the  precise  physical  nature  of  the 
incremental  stresses  S..  and.  S...  Since  S-.  represents  the  incremental 

o  i j  t  ij  o  IJ 

2  a  Piola-Kirchhoff  stress  tensor  at  time  t  referred  to  the  configuration 

at  time  o,  it  is  a  convected  quantity  deforming  together  with  the  material 

element  of  a  body.  The  stress  quantity  .  S..  has  similar  definition  as  S.  • 

t  ij  o  ij 

except  that  it  is  referred  to  the  updated  configuration  at  time  t.  That 
is  to  say,  is  not  deforming  with  the  material  fibers.  For  an  ortho¬ 

tropic  or  anisotropic  material  undergoing  large  displacements,  the  use  of 
t  S . j  in  the  constitutive  description  will  not  be  appropriate  unless  the 
effect  of  material  axes  rotations  is  accounted  for.  Such  option  is  not 
available  in  ADINA.  Nevertheless,  the  program  does  impose  the  necessary 
restriction  that  only  the  TL  formulation  is  permitted  for  the  case  of 
large  deformation  analysis  of  orthotropic  materials. 

Another  important  note  is  the  limitation  of  the  program  in  dealing 
with  material  incompressibility.  In  engineering  application,  the  elastic 
material  most  frequently  used  for  large  strain  analysis  is  the  Mooney- 
Rivlin  (or  rubber-like)  material.  In  order  to  handle  the  general  case 
nonlinear  analysis  of  such  material,  an  incompressible  finite  element 
formulation  [  7]  is  necessary.  Moreover,  in  large  strain  elastic-plastic 
analysis  such  as  metal-forming  problems,  plastic  incompressibility  condition 
introduces  two  numerical  difficulties  in  the  finite  element  analysis: 
i)  the  element  becomes  overly  constrained  causing  unnecessary  stiffening 
effect  [31]  ,  and  ii)  small  time  (or  loading)  increment  is  necessary 
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in  order  to  obtain  convergent  solution  [29]  .  From  this  discussion,  it 
appears  that  ADINA  will  have  greater  versatility  if  a  finite  element 
formulation  with  incompressibility  constraint  is  included  in  the  future 
extension  of  the  code. 

As  mentioned  earlier,  the  material  models  in  ADINA  are  element- 
dependent,  both  the  2/D  and  3/D  continuum  elements  have  a  fairly  extensive 
list  of  material  models.  From  the  theoretical  consideration,  the  material 
models  included  in  ADINA  can  be  classified  in  two  categories:  i)  classical 
models,  which  conform  to  the  basic  principles  in  mechanics,  and  ii)  empirical 
models,  which  are  largely  based  on  experimental  observations  and/or  con¬ 
jectures.  In  accordance  with  these  classifications,  the  material  models 
are  listed  as  the  following: 

Classical  Models 

Isotropic  linear  elastic 
Orthotropic  linear  elastic 
Isotropic  thermo-elastic 
Mooney-Rivl in  material 

Drucker-Prager  model  (Elastic-perfectly  plastic  material) 
von  Mises  elastic-plastic  with  isotropic  or 
kinematic  hardening 
Thermo-elastic-plastic  and  creep 

Empirical  Models 


Curved  description  model 
Concrete  model 
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In  addition  to  the  above  discussion,  some  further  comments  can  be 

made: 

1)  The  von  Mises  elastic-plastic  model  with  either  isotropic  or 
kinematic  hardening  is  restricted  to  materials  with  bilinear 
stress-strain  curve  in  uniaxial  response.  This  limitation  is 
unnecessary  for  the  isotropic  hardening  materials  and  the  code 
could  be  modified  to  include  the  nonlinear  stress-strain  response. 

2)  The  curve  description  model  was  intended  for  simulating 
the  nonlinear  response  of  geological  materials  [ 32  ]  .  In 
formulating  the  constitutive  relations,  a  number  of  assumptions 
were  made.  First,  the  incremental  stresses  and  strains  are 
related  by  instantaneous  moduli  in  shear  and  bulk,  respectively. 


where 


t 


G 


2  ‘G  5*3 

3  tK  e 

m 

Incremental  deviatoric  stresses 
Incremental  deviatoric  strains 
Incremental  mean  stress 
Incremental  mean  strain 

Instantaneous  bulk  modulus,  which  is  a  function 
of  loading  condition  and  volume  strain.  One 
set  of  input  data  is  required  for  loading  and 
another  set  for  unloading. 

Instantaneous  shear  modulus,  which  is  also  a 
function  of  loading  condition  and  volume  strain. 
Only  one  set  of  input  data  is  necessary  for  loading 
condition.  The  values  for  unloading  are  scaled  by 
the  stress-strain  data  in  bulk. 


(25) 

(26) 
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The  loading  and  unloading  conditions  of  the  material  are  specified 

from  a  strain  criterion:  If  the  calculated  te  <  e  .  (a  minimum  mean 

m  —  min 

strain),  the  material  is  under  loading,  and  if  ^e  >  e  .  ,  the  material 

m  mi  n 

is  under  unloading.  In  addition,  the  material  may  weaken  or  crack  when 
the  maximum  principal  stresses  exceed  certain  limit. 

It  is  important  for  the  user  to  note  three  possible  limitations 
of  this  model: 

i)  The  material  behavior  is  uniquely  determined  by  the 
magnitude  of  volume  strain  without  regard  to  its  deformation 
history  although  different  stress-strain  curves  are  used  for 
loading  and  unloading. 

ii)  The  consti  tutive  relationship  assumes  that  the  principal 
axes  of  the  stress  increment  coincide  with  those  of  the  strain 
increment.  This  is  generally  not  the  case  due  to  material  - 
induced  anisotropy. 

iii)  The  relationships  in  Eqs.  (25)  and  (26)  are  equivalent  to 
a  non-associated  flow  rule  which  is  known  for  lack  of  solution 
uniqueness,  excepted  in  a  very  specialized  case  [  33  1 . 

3)  The  concrete  model  in  ADINA,  although  involves  some  sophisticated 
assumptions  [34],  is  very  similar  to  the  curve  description  model.  There¬ 
fore,  this  model  is  also  subject  to  the  limitations  cited  above.  Further¬ 
more,  to  model  the  nonl inear-fracturing  behavior  of  concrete  materials, 
a  number  of  different  constitutive  theories  have  been  proposed  in  the 
literature  [35-39]  and  the  user  should  make  proper  judgement  on  the 
selection  of  the  models  that  are  available. 
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3.  Finite  Element  Idealization 

In  ADINA,  the  truss,  2/D  continuum  and  3/D  continuum  elements 
are  the  isoparametric  family  except  the  beam  element  for  which  inde¬ 
pendent  interpolation  functions  are  used  to  approximate  the  displace¬ 
ment  field.  The  use  of  isoparametric  elements  has  several  apparent 
advantages  [7,40]: 

i)  The  elements  satisfy  the  convergence  requirements  of  finite 
elements.  That  is,  rigid  body  modes,  constant  strain  behavior 
and  inter-element  compatibility  are  preserved. 

ii)  Various  structural  configurations  and  curved  boundaries  can 
be  modeled. 

iii)  Different  material  models  can  easily  be  incorporated  into 
the  stiffness  matrix  by  numerical  integration. 

iv)  The  elements  give  better  stress  field  as  compared  to  the 
triangular  or  tetrahedron  elements. 

For  a  quick  reference,  the  element  types  vs.  degrees  of  freedom 
and  order  of  numerical  integration  are  shown  in  Table  4.  A  brief  out¬ 
line  on  the  finite  element  approximations  for  each  element  type  is 
given  below. 

3.1  Truss  Element 

It  represents  a  structural  member  capable  of  transmitting  the 
axial  force  only.  The  element  may  have  2,  3,  or  4  nodes,  and  each 
node  has  3-degrees  of  freedom  in  translation  (Fig.  1).  The  user  may 
also  have  the  option  to  use  the  truss  element  for  axisymmetric  analysis 
in  which  case  only  one  node  is  defined  in  the  YZ-plane.  The  inter- 
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Table  4.  Order  of  Numerical  Integration 


Element  Type 

No.  of  Nodes 

per  Element 

DOF  per 

Node 

Order  of  Numerical 

Integration 

Truss 

2  -  4 

3 

2  -  4 

2/D  Continuum 

4-8 

2 

2-4 

3/D  Continuum 

8  -  21 

3 

2  -  4 

Beam 


2 


6 


r-direction:  1-7 
s-direction:  1-7 
t-direction:  1-8 


polation  functions  used  for  the  element  formulation  can  be  found  from 
references  [41,42].  Since  the  cross  sectional  area  was  assumed  to 
remain  constant  during  deformation,  the  element  can  be  applied  only 
to  large  displacement  analysis  with  small  strains.  Large  deformation 
analysis  is  based  on  the  UL  formulation. 

3.2  2/D  Continuum  Element 

It  is  an  isoparametric  element  with  4-8  nodes,  and  2-degrees 
of  freedom  at  each  node  as  seen  in  Fig.  2.  The  interpolation  func¬ 
tions  of  this  element  can  be  found  from  references  [7,41],  Although 
the  user  has  the  option  to  choose  any  number  of  nodes  between  4  and 
8  to  describe  an  element,  use  of  5-  node  and  7-  node  elements  should 
be  avoided  due  to  the  incompatible  displacement  field  in  the  element 
which  may  cause  numerical  difficulty  [43]. 

The  element  stiffness,  consistent  mass  and  internal  force  are 
evaluated  by  numerical  integration  with  Gauss  quadrature  formulae. 

The  user  must  specify  the  order  of  numerical  integration  varying  frqm 
2  to  4.  Generally  for  linear  elastic  material,  an  integration  order 
of  2  is  sufficient.  However,  for  nonlinear  material  in  which  high 
stress  gradients  occur,  an  integration  order  of  3  is  recommended.  The 
program  sets  the  integration  order  of  3  for  the  evaluation  of  consis¬ 
tent  mass  matrix. 

A  triangular  element  can  be  obtained  by  degenerating  a  4-  node 
or  8-  node  quadrilateral  element  as  shown  in  Fig. 3.  The  use  of  degener¬ 
ated  8-node  triangular  element  should  be  avoided  since  the  element  stiff- 


ness  is  biased  depending  on  the  manner  that  the  local  nodes  are  com¬ 
bined  [7,40].  However,  the  degenerated  8-  node  element  with  quarter 
point  has  proven  to  be  useful  to  model  the  crack  tip  singularity  in 
fracture  mechanics  [44]. 

3.3  3/D  Continuum  Element 

The  numerical  approximation  of  this  element  [7,41]  is  basically 
the  same  as  the  2/D  element,  therefore,  most  of  the  discussion  given 
in  the  above  is  also  applicable  here.  The  element  consists  of  8  - 
20  nodes  with  3-  degrees  of  freedom  at  each  node  (Fig.  4).  Gauss 
quadriture  formulae  is  used  in  numerical  integration  with  the  integra¬ 
tion  order  varying  between  2  and  4.  The  user  may  specify  the  inte¬ 
gration  order  in  the  surface  and  over  the  thickness  of  the  element  inde¬ 
pendently.  Since  the  computational  cost  for  evaluating  the  stiffness 
matrix  of  a  3/D  element  is  rather  high  and  it  is  closely  associated 
with  the  number  of  integration  points  chosen,  the  user  must  pay  extra 
care  in  making  such  a  decision.  In  most  cases,  an  integration  order 
of  2  is  sufficiently  accurate.  However,  in  case  of  high  stress  gradi¬ 
ents  or  nonlinear  materials,  an  integration  order  of  3  is  necessary. 

If  this  element  is  to  be  applied  to  the  bending  analysis  of  a  thin 
plate  or  a  thin  shell,  reduced  order  of  numerical  integration  must 
be  used  in  the  plane  of  the  element  in  order  to  obtain  any  meaningful 
results.  This  is  further  discussed  in  Part  II.  in  the  report. 
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3.4  Beam  Element 

Differing  from  the  previous  three  elements,  the  beam  element  is 
not  a  member  of  the  isoparametric  family.  The  element  is  described 
by  two  nodes  at  the  ends  with  6-  degrees  of  freedom  at  each  node: 

3-  translations  and  3-  rotations  (Fig.  5).  In  the  finite  element  formu¬ 
lation,  the  following  assumptions  were  made: 

i)  Only  prismatic  straight  beam  was  considered, 

ii)  Plane  sections  of  the  beam  remain  plane  during  deformation, 
but  not  necessarily  perpendicular  to  the  neutral  axis,  i.e., 
a  constant  shear  is  allowed. 

iii)  The  beam  undergoes  large  rotations  but  small  strains, 

i v )  A  cubic  variation  in  bending  displacement  and  linear  variation 

in  axial  and  torsional  deformations, 

v)  The  UL  formulation  was  adopted. 

The  interpolation  functions  for  the  displacement  field  are  based 
on  those  in  reference  [42]  and  the  derivation  of  nonlinear  stiffness 
matrices  are  given  in  [45].  Numerical  integration  with  Newton-Cotes 
formulas  is  used  for  the  evaluation  of  nonlinear  element  stiffness 
matrix  and  internal  force  vector.  The  integration  order  varies  from 
5  to  7  in  the  longitudinal  direction  of  the  beam  and  3  to  7  in  the 
transverse  directions.  However,  for  a  linear  elastic  beam  under  small 
deformation,  the  stiffness  matrix  is  evaluated  by  close  form  formu¬ 
lations  [42,46]. 
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4.  Static  Solution  Method 

The  solution  method  employed  in  ADINA  for  solving  nonlinear  static 
problems  is  based  on  the  incremental  loading  with  tangent  stiffness 
approach  and  equilibrium  iterations,  which  corresponds  to  the  Modified 
Newton-Raphson  (MNR)  procedure[7, 23, 47-49] .  The  user  has  the  option 
to  request  the  stiffness  reformation  at  every  loading  increment  or 
at  specified  loading  increments.  The  user  also  controls  the  number 
of  equilibrium  iterations  for  convergence.  If  this  number  is  not 
specified,  the  program  set  it  to  be  "15"  internally.  The  apparent 
advantage  of  this  method  is  that  it  is  computationally  efficient  for 
solving  many  practical  nonlinear  problems.  Fig.  6  shows  the  converg¬ 
ence  characteristics  of  the  MNR  method  in  the  load-deflection  plot  of 
a  structure  exhibiting  softening  behavior.  It  is  seen  that  if  the  load 
increment  is  sufficiently  small,  solution  will  converge  monotonically. 
Convergence  of  the  solution  is  determined  in  the  program  by  checking 
the  calculated  displacements  during  iterations,  i.e. 

<  RTOL  (27) 

where  ||  ||  denotes  the  Euclidean  norm  and 

AU^  =  Incremental  displacement  calculated  up  to  the  i-th 
iteration  cycle. 

y(t+i-l)  _  displacement  calculated  up  to  and  including  the 

current  iteration  cycle. 

RTOL  =  A  convergence  tolerance  specified  by  input,  usually 
in  the  order  of  0.01  or  smaller. 
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In  addition,  the  solution  from  iterations  is  declared  to  diverge  if 
the  norm  of  the  out-of-balance  load  is  larger  than  the  norm  of  the 
incremental  load  for  the  current  solution  step. 

Generally  the  MNR  method  does  not  give  convergent  solution  under 
three  conditions: 

1)  Load  increment  is  too  large. 

2)  Problems  involve  elastic-plastic  deformations  when  elastic 
unloading  occurs.  One  way  to  circumvene  Shis  problem  is 
to  either  use  elastic  material  law  for  the  elements  under¬ 
going  unloading  or  specify  a  very  small  load  increment  and 
suspend  the  equilibrium  iteration  just  to  overcome  the 
elastic  unloading. 

3)  The  structure  is  becoming  stiffened  during  deformation.  In 
this  case,  the  solution  may  or  may  not  converge  depending 
on  the  load  increment  specified  or  the  load-deformation 
characteristics  as  seen  in  Fig.  7-a  and  7-b.  Even  if  the 
convergence  does  occur,  the  iterative  displacements  and 
residual  forces  are  oscillatory  and  the  convergence  rate 

is  generally  slow. 

Further  discussion  on  the  numerical  characteristics  of  the  MNR 
method  will  be  given  in  Part II  of  this  report. 


...  .Ci _ 
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5.  Dynamic  Analysis  Method 

With  the  inclusion  of  the  inertia  effect  in  the  virtual  work  equa¬ 
tion,  i.e.  Eq.  (11)  or  (12),  and  after  the  finite  element  discretiza¬ 
tion,  the  equations  of  motion  of  a  structure  corresponding  to  either 
TL  formulation  or  the  UL  formulation  can  be  written  in  a  typical  matrix 
form  at  t  +  r 


Mt+Tii  +  Ct+TU  +  Kt+TU  =  t+TR  - 

where 


M 

C 

K 

U 

U 

U 

t+TR 


Mass  matrix 

Damping  matrix 

Structural  stiffness  matrix 

Incremental  nodal  displacement  vector 

Incremental  nodal  velocity  vector 

Incremental  nodal  acceleration  vector 

Total  external  nodal  force  vector  at  time  t+T 


(28) 


i 

* 


r 

.4 


Total  internal  nodal  force  vector  at  time  t. 


For  nonlinear  problems,  the  matrix  K  represents  the  tangent  stiffness 
matrix  of  the  finite  element  idealization  which  may  include  the 
effects  due  to  geometric  and  material  nonlinearities.  The  mass 
and  damping  matrices  in  the  above  equations  are  the  same  as  those 
of  a  linear  structure. 

In  ADINA,  Eq.  (28)  is  solved  by  using  the  direct  integration 
method  and  both  the  implicit  and  explicit  integration  schemes  with 
the  following  operators  are  available  in  the  program: 

Implicit  Integration 

i)  Newmark  8  -  method  [50] 
ii)  Wilson  6  -  method  [51] 

Explicit  Integration 

i)  Central  -  difference  method  [52]. 

The  stability  and  accuracy  of  the  aforementioned  operators  have 

been  investigated  quite  extensively  for  linear  structural  dynamics  by 

a  number  of  investigators  [53-60].  These  studies  showed  that  both  the 

Wilson  9-method  and  Newmark  8  -  method  (y  =  H,  8  =  h)  are  unconditionally 

stable,  whereas  the  central  difference  method  is  only  conditionally 

stable.  The  stability  condition  for  the  central  difference  method  is 

that  the  time  step  At  must  be  smaller  than  a  critical  value  At  [41]  i.e., 

cr 
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where 

Tn  =  smallest  vibrational  period  of  the  finite  element 
system  to  be  considered. 

From  the  study  made  by  Goudreau  and  Taylor  [57],  the  Wilson  method  was 
found  to  introduce  considerable  damping  to  the  structure  and  the  Newmark 
average  acceleration  scheme  does  not.  They  concluded  that  the  use  of 
Newmark  method  requires  less  than  half  as  many  time  steps  as  the 
Wilson  method. 

While  most  experience  on  the  stability  and  solution  accuracy  of 
the  approximate  operators  was  gained  for  the  linear  dynamic  problems, 
little  is  known  on  the  application  of  these  operators  to  nonlinear 
problems.  Weeks  [61]  has  evaluated  the  Houbolt  and  Newmark  operators 
for  a  one- degree-of- freedom  geometrically  nonlinear  problem  and  found 
that  both  operators  gave  a  satisfactory  solution  when  the  full  Newton-Raphson 
iteration  scheme  is  used.  This  result  was  also  substantiated  by  McNamara 
[62]  from  his  numerical  investigations.  From  the  energy  considerations, 
the  Newmark  method  with  8  _>  k  was  also  found  to  be  unconditionally  stable 
when  applied  to  nonlinear  structural  dynamics  [63-65].  Additional 
computational  aspects  in  regard  to  the  application  of  the  direct  inte¬ 
gration  method  are  discussed  in  [66-68]. 

Despite  all  the  progress  that  has  :>een  made  in  structural  dynamics, 
application  of  ADINA  or  any  other  computer  codes  to  the  analysis 
of  nonlinear  sturctural  problems  is  a  difficult  task.  The  analysts 
are  always  confronted  with  two  major  issues: 

i)  A  time  step  which  can  give  convergent  and  reliable  solutions 
ii)  Computational  costs 
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Unfortunately,  there  are,  as  yet,  no  simple  answers  to  these  issues. 

From  the  literature,  some  of  the  points  which  might  be  helpful  to  the 
analysts  are  summarized  as  the  following: 

1)  The  apparent  advantage  of  the  explicit  integration  scheme  over  the 
implicit  scheme  is  that  for  a  lumped  mass  system  it  does  not  require 
the  assemblage  of  structural  and  mass  matrices,  and  matrix  inversion 
which  represent  a  major  computational  cost  in  the  finite  element  analysis. 
However,  this  aspect  is  hampered  by  the  fact  that  the  explicit  scheme 

is  only  conditionally  stable  in  numerical  integration  and  hence  it 
requires  smaller  time  steps  as  compared  to  the  use  of  implicit  schemes. 

2)  Some  general  discussions  on  the  choice  of  the  explicit  and  implicit 
schemes  were  given,  for  example,  in  [54,  56,  66,68].  For  wave  propagation 
problems  where  high  frequency  response  of  the  structure  is  more  important 
and  the  discontinuity  in  velocity  or  acceleration  persists,  the  explicit 
integration  scheme  will  be  the  obvious  choice.  On  the  other  hand,  for 
structural  dynamic  problems  where  the  response  is  predominated  by  the 
fundamental  frequencies  of  the  system,  the  implicit  scheme  is  preferred. 

3)  When  comparing  the  Newmark  8-method  vs.  Wilson  6-method,  in 
most  cases  both  methods  give  the  same  results  and  have  the  same  compu¬ 
tational  characteristics  [41].  However,  the  Newmark  method  conserves 
the  total  energy  for  an  undamped  system  whereas  the  Wilson  method  does 
not  [57].  Therefore,  the  Newmark  method  seems  to  be  superior.  Further¬ 
more,  this  method  was  proven  to  be  unconditionally  stable  when  applied 
to  nonlinear  dynamic  problems. 


4)  For  implicit  integration,  the  step  At  should  be  smaller  than 
(1/20;!^,  where  Tn  is  the  maximum  period  of  the  nonlinear  system  [69]. 

5)  Computationally,  it  is  more  efficient  to  invoke  equilibrium 
iterations  in  every  loading  step  while  the  stiffness  is  updated 
at^every  2- 3.. losing  step. 
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IV.  PROGRAMING  DESCRIPTION 

In  its  early  inception,  a  general  purpose  finite  element  struc¬ 
tural  software  mostly  represents  an  extension  of  special  purpose 
finite  element  codes  by  utilizing  certain  common  features: 

1)  Reading  and  processing  of  input  data 

2)  Assemblage  of  the  master  matrix 

3)  Solution  of  large  systems  of  matrix  equations 

4)  Performing  certain  analysis  procedures  (i.e.  static  or 
dynamic,  linear  or  nonlinear  analysis) 

5)  Printing  of  output,  etc. 

With  these  features  included,  various  element  types  and  material  de¬ 
scriptions  can  be  added  into  the  software  to  form  a  general  purpose 
code  for  conducting  a  wide  range  of  structural  analysis.  However, 
a  large  scale  program  without  proper  coding  organization  will  soon 
become  obsolete  as  it  is  expanded  significantly  beyond  its  original 
capability.  More  specifically,  a  poorly  designed  general  purpose 
software  may  face  the  following  potential  problems: 

1)  It  takes  much  longer  computer  time  due  to  many  conditional 
statements  and  unnecessary  data  transfers  between  the  fast 
core  and  low  speed  storages. 

2)  It  requires  excessive  computer  storage  for  poor  data 
management. 

3)  Coding  changes  become  difficult  without  interfering  with 
other  parts  of  the  software. 


To  alleviate  some  of  these  problems,  a  large  scale  general  purpose 
software  must  be  written  on  a  modular  basis  and  only  a  few  modules 
are  involved  for  a  particular  analysis.  Moreover,  all  modules  should 
be  organized  in  a  manner  that  they  can  be  modified  independently  or 
changed  without  interfering  with  each  other.  In  view  of  these  con¬ 
siderations,  the  modern  design  of  general  purpose  finite  element  soft¬ 
ware  has  evolved  into  a  characteristic  form.  That  is,  the  entire  soft¬ 
ware  may  be  divided  into  three  major  parts: 

1)  Input  phase  (pre-processor) 

2)  Analysis  phase  (main  program) 

3)  Output  phase  (post-processor) 

Since  in  ADINA,  the  output  phase  is  a  part  of  the  analysis  phase  in 
terms  of  its  coding  organization,  both  items  2  and  3  will  be  reviewed 
as  one  unit.  To  provide  a  basis  for  discussion  and  to  assist  the 
ADINA  user  to  follow  its  macro  command,  a  flow  diagram  for  the  main 
module  "AAMAIN"  is  included  in  the  Appendix  Al. 

Also  included  in  this  chapter  is  a  dis;ussion  on  the  use  of  con¬ 
trol  variables  and  data  structure  adopted  in  the  code. 

1.  Input  Phase 

The  major  function  of  this  phase  is  to  read  and/or  generate 
input  data  which  are  necessary  to  define  a  structural  problem.  These 
i ncl ude : 

1)  Read  the  definitions  of  master  control  variables, 

time  history  and  integration  method,  printing  options,  etc. 
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Table  5. 

Input  Sequence 

Sequence 

Subroutines 

Flaq 

Input  Description 

1 

ADINI 

IND  =  0 

Definition  of  master  control 
variables  for  solution  options, 
time-history,  and  nodal  print 
data. 

2 

ELCAL 

IND  =  0 

Material  data  and  element 
definitions. 

3 

ADINI 

IND  =  2 

Applied  loads:  Concentrated 
forces  or  pressure  load. 
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2)  Read  nodal  constraint  condition  and  coordinates,  then 
generate  the  equation  number  matrix  (ID-array)  and  store 
it  on  the  low  speed  device,  unit  #8. 

3)  Read  the  initial  conditions  and  initialize  the  nodal 
displacements,  velocities  and  accelerations.  Write  this 
information  on  unit  #8. 

4)  Read  and  generate  element  data:  material  properties 
and  element-nodal  number  relations.  From  these,  form  the 
element  connection  arrays  (LM)  and  store  on:  unit  #1  -  for 
linear  elements,  unit  #2  -  for  nonlinear  elements. 

The  above  executions  are  done  in  three  stages  by  calling  the  sub¬ 
routines  ADINIand  ELCAL.  The  calling  sequence  is  controlled  in 
the  main  driver  AAMAIN  (as  shown  in  Al)  by  use  of  a  control  vari¬ 
able  IND  (=0  or  2).  This  interrelationship  is  shown  in  Table  5. 
For  further  references,  the  subprograms  called  by  AD IN I  and  ELCAL 
are  shown  in  sections  A2  and  A3  in  the  Appendix  A,  respectively. 

It  is  noted  that  all  input  data  must  be  provided  according 
to  the  fixed  format  strictly  following  the  sequential  order  in 
the  user's  manual.  Reading  and  calculation  of  load  data  are  made 
after  the  assemblage  of  linear  structural  stiffness  (see  flow 
diagram  Al).  This  calling  sequence  cannot  be  easily  reversed 
due  to  the  dynamic  allocations  of  array  storage  for  the  load 
calculations  coded  in  the  subroutine  AD I N I . 


2.  Analysis  Phase 


The  analysis  phase  represents  the  main  package  of  the  Pro¬ 
gram  which  carries  out  the  actual  problem  solving.  The  composition 
of  this  phase  is  generally  made  up  by  three  libraries: 

1)  Procedure  library  -  which  consists  of  the  solution 
options  available  to  the  user. 

2)  Geometric  library  -  which  includes  various  element 
types  for  the  geometric  modeling  of  a  structure. 

3)  Material  library  -  which  contains  various  constitutive 
models  for  simulating  the  material  behavior. 

Each  of  the  above  categories  is  reviewed  as  below. 

2.1  Procedure  Library 

The  analysis  procedures  are  primarily  controlled  by  the  main 
driver  "AAMAIN".  The  options  available  include: 

1)  Linear  static  analysis 

2)  Linear  dynamic  analysis  (direct  step-by-step  integration) 

3)  Nonlinear  static  analysis 

4)  Nonlinear  dynamic  analysis 

5)  Frequency  analysis 

6)  For  dynamic  analysis,  the  choice  of  integration  operators 
Wilson  e,  Newmark  or  central  difference  method. 

7)  In  core  or  out-of-core  solution  (determined  internally). 

8)  Reformation  of  structural  stiffness 

9)  Iteration  option  for  nonlinear  static  or  dynamic  analysis 


WiialtiCli 
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Decision  on  the  analysis  options  is  made  in  AAMAIN  by  checking 
the  pre-a ssigned  values  of  the  control  variables  (or  flags)  and 
subsequently,  various  analysis  modules  are  called  to  carry  out 
the  analysis  options  requested  by  the  user.  An  overview  of  the 
analysis  procedures  is  shown  in  Fig.  8  and  a  detailed  flow  diagram 
for  AAMAIN  is  given  in  Section  A1  of  the  Appendix  A. 

As  seen  from  Fig.  8  and  Section  A1 ,  the  modules  corresponding 
to  various  analysis  procedures  were  arranged  in  a  serial  form 
in  AAMAIN.  For  any  desired  analysis  option,  a  specific  path  is 
channeled  through  by  combining  appropriate  modules  or  subroutines 
with  the  use  of  conditional  statements.  For  discussion,  the  macro 
comnand  of  the  analysis  procedures  (1)  -  (4)  are  listed  in  Tables  6-9 
respectively.  The  calling  sequences  o'  the  modules  for  dif¬ 
ferent  analysis  procedures  are  quite  similar  with  some  minor  vari¬ 
ations.  For  example,  to  compare  the  algorithm  difference  between 
the  linear  static  analysis  and  linear  dynamic  analysis,  the  macro 
command  is  almost  identical  except  that  three  modules  perform 
somewhat  different  calculations  as  seen  in  the  following: 


Module 

Linear  Static  Analysis 

Linear  Dynamic  Analysis 

ASSEM 

Form  global  stiffness 

Form  effective  global 
stiffness,  mass  and 

damping  matrices. 

LOADEF 

No  calculations  are 

Form  effective  load  vec¬ 

performed 

tors  due  to  dynamics 

NEWDAV 

- 

Update  nodal  displacements, 
velocities  and  accelerations 

FORM  LINEAR  STRUCTURAL 
STIFFNESS 
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Table  6.  Macro  Command  for  Linear  Static  Analysis 


Module 


Description 


ADDRES 

SBLOCK 

ASSEM 

ADINI 

WRITE 

RSTART 


RSTART 


Determine  the  locations  of  diagonal  elements 
of  global  stiffness 

Determine  block  number  and  size  for  in-core 
or  out-of-core  solution 

Calculate  element  stiffness  and  assemble 
global  stiffness 

Construct  load  vectors 

Print  out  initial  displacements 

For  a  restart  job,  read  displacement  vector 
from  tape  #8  (optional) 

Triangularize  linear  stiffness  matrix 

Obtain  solution  by  load  reduction  and  back 
substitution 

Update  nodal  displacements 

Print  out  nodal  displacements 

Calculate  element  strains  and  stresses  at 
Gauss  points 

Write  displacement  vector  on  tape  #8  for  a 
restart  job  (optional) 


Time  Looi 
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Table  7  Macro  Command  for  Linear  Dynamic  Analysis 
and  Frequency  Analysis 


Module 


Description 


ADDRES 

(1) 

SBLOCK 

(1) 

ASSEM  Calculate  element  stiffness  and  mass  matrices 

Form  effective  global  stiffness,  mass  and 
damping  matrices. 

ADINI  (1) 

WRITE  Print  out  initial  displacements,  velocities  and 

accelerations 


RSTART 

ASSEM 

C0LS0L 

LOADEF 

C0LS0L 

NEWDAV 

WRITE 

STRESS 


For  a  restart  job,  read  displacements,  velocities 
and  accelerations  from  tape  #8  (optional) 

Calculate  frequencies  and  mode  shapes  (2) 

(1) 

Form  effective  load  vectors  due  to  dynamic  effects 

(1) 

Update  nodal  displacements,  velocities  and 
accelerations 

Print  out  nodal  displacements,  velocities  and 
accelerations 


(1) 


RSTART 


Write  displacements,  velocities  and  accelerations 

on  tape  #8  for  a  restart  job  (optional)  ir 


XT)  Same  as  Table  7 

(2)  For  frequency  analysis  option 


Load  Increment  Loop 


m 
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Table  8 

Macro  Command  for  Nonlinear  Static  Analysis 

Module 

Description 

ADORES 

(1) 

SBLOCK 

(1) 

ASSEM 

Perform  the  same  calculations  as  in  Table  7 
for  linear  element  group 

ADINI 

(1) 

WRITE 

(1) 

RSTART 

(1) 

ASSEM 

Calculate  tangent  element  stiffness  and 
internal  nodal  forces. 

Assemble  global  stiffness  and  form 
incremental  load 

COLSOL 

Obtain  incremental  displacements 

EQUIT 

Perform  equilibrium  iterations 

NEWDAV 

(1) 

WRITE 

(1) 

STRESS 

(1) 

RSTART 

(1) 

JJ)  Same  as  Table  7 


s4 

o 


i 
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Table  9. 

Macro  Command  for  Nonlinear  Dynamic  Analysis 

Module 

Description 

ADORES 

(1) 

SBLOCK 

(1) 

ASSEM 

Form  mass  and  damping  matrices 

Form  stiffness  m<trix  for  linear  element 

ADINI 

(1) 

WRITE 

(2) 

RSTART 

(2) 

LOADEF 

(2) 

ASSEM 

Calculate  tangent  element  stiffness  and 
internal  nodal  force 

Assemble  global  stiffness  and  incremental 

COLSOL 

(3) 

EQUIT 

(3) 

NEWDAV 

(2) 

WRITE 

(2) 

•STRESS 

(1) 
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It  is  noted  that  for  the  static  analysis  the  module  LOADEF  is 
still  entered  and  then  immediately  returned  without  performing 
any  calculations.  Another  point  of  difference  is  that  the  dyna¬ 
mic  analysis  is  repeated  over  a  time  loop  as  indicated  in  Table  7 
whereas  the  static  analysis  may  loop  over  the  multiple  load  steps. 

By  comparing  Tables  8  and  9  with  Tables  6  and  7,  the  major  dif¬ 
ferences  between  the  nonlinear  and  linear  analyses  (for  both 
static  and  dynamic)  are: 

1)  Formation  of  tangent  stiffness  in  ASSEM 

2)  Performing  equilibrium  iterations  in  EQUIT 

In  conjunction  with  the  aforementioned  macro  command,  the  central 
memory  for  data  storage  and  work  space  is  utilized  by  defining  a 
master  array  in  the  blank  common  block.  This  master  array  is  further 
divided  sequentially  into  subregions  for  arrays  which  have  variable 
lengths  (also  called  dynamic  allocations)  depending  on  the  storage 
requirements  of  the  modules  involved.  In  the  main  driver  AAMAIN,  array 
sizes  are  defined  prior  to  the  analysis  modules  that  are  called 
and  hence  the  pointer  for  dynamic  allocations  vary  from  one  stage 
of  the  analysis  to  the  other.  Because  of  this  relationships, 
we  may  draw  the  following  characteristic  statements: 

1)  The  procedure  library  is  commanded  in  AAMAIN  by  use 
of  control  variables. 

2)  Array  allocations  in  the  blank  common  block  vary  with 
the  algorithm  procedures  and  the  manner  of  partitioning 
the  master  array  for  various  analysis  steps  are  inter¬ 
related. 
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3)  Coding  organization  of  the  procedure  library  does  not 
conform  to  the  modular  concept.  Consequently,  it  is 
rather  difficult  to  alter  this  part  of  the  coding  in 
ADINA  without  upsetting  the  existing  analysis  procedures. 
For  completeness,  the  flow  diagrams  for  the  following  compu¬ 
tations  are  also  included  in  the  Appendix  A. 

A4  -  Frequency  Analysis 
A5  -  Calculation  of  Effective  Load 
A6  -  Assembly  of  Structural  Stiffness 
A7  -  Equilibrium  Iteration 
A8  -  Stress  Calculation 

2.2  Geometric  Library 

The  Geometric  library  is  composed  of  various  element  type 
which  are  used  to  model  the  geometric  configuration  and  response 
of  a  structure  under  given  loads.  This  library  is  commanded  by 
the  routine  "ELEMNT"  with  ELCAL  or  ASSEM  being  a  connecting  routine 
between  ELEMNT  and  the  main  driver  A AM A IN  as  shown  in  Fig.  9. 

The  routine  ELEMNT  then  calls  for  four  element  control  rountines, 
i  .e. 

1)  TRUSS  -  Truss  element 

2)  TODMFE  -  2/D  Continuum  element 

3)  THREDM  -  3/D  Continuum  element 

4)  BEAM 


Beam  element 
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*  ELCAL  -  for  input  of  element  data. 

ASSEM  -  for  element  stiffness  or  load  calculation. 


Fig.  9.  Element  Library 
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The  flow  diagrams  for  the  above  element  packages  are  further  illus-  j  : 

trated  in  sections  A9  -  A12  of  the  Appendix  A.  Each  of  the  above  j 

routines  dynamically  allocates  the  in-core  storage  space  for  data 

processing  of  that  element.  Then  the  element  control  routine 

calls  for  the  lower  level  routines  tc  perform  the  following  functions: 

i  < 

1)  Read  and  generate  element  definitions  and  material  > 

properties.  1 

2)  Form  the  linear  element  stiffness  matrix. 

3)  In  case  of  dynamic  analysis,  form  the  lumped  (or  consis¬ 
tent)  mass  matrix.  I 

4)  Form  the  nonlinear  (material  or/and  geometric)  stiffness  , 

matrix  and  internal  nodal  force  vector. 

5)  Calculate  element  strains  and  stresses  and  print  out  j  j 

these  quantities. 

All  the  above  functions  are  coded  in  the  sub-module,  e.g.  TDFE 
for  2/D  continuum  element.  The  calling  sequence  of  the  element 
module  to  perform  the  aforementioned  function  is  shown  in  Table  10. 

Execution  decision  on  these  functions  are  made  by  checking  the 

J 

present  values  of  two  control  variables: 

IND  -  Flag  for  analysis  phase 

KPRI  -  Flag  for  strain  and  stress  calculations  ^  I 

Among  the  four  element  types  include!  in  ADINA,  the  coding  for 
the  2/D  element  (TDFE)  and  3/D  el  erne  it  (THDFE)  was  written  in 
a  parallel  pattern.  However,  the  coding  for  the  truss  and  beam 
elements  does  not  appear  to  be  as  well  organized  as  the  continuum  j1 

elements.  h 


■nr 
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Table  10.  Calling  Sequence  of  the  Element  Module 


Analysis  Stage 

Callinq  Routines 

Description 

Input 

ELCAL 

Read  element  data  and 
material  properties 

Formation  of  Element 
Matrices 

ASSEM 

Formation  of  element 
stiffness  and  mass 
matrices,  and  calculate 
internal  nodal  force 
vector. 

Equilibrium  Iteration 

EQUIT 

Calculate  internal  nodal 
forces  for  nonlinear 
elements. 

Stress  Calculations 

STRESS 

Calculate  element  strains 
and  stresses 

One  important  characteristics  of  the  element  module  is  that 
each  element  type  is  a  self-contained  unit  which  is  completely 
independent  from  other  elements.  Therefore,  modification  or  addi¬ 
tion  of  an  element  to  the  program  is  a  relatively  simple  task. 

2.3  Material  Library 

As  noted  previously,  the  material  library  in  ADINA  was 
coded  as  a  subsidiary  of  the  element  module.  That  is,  each  element 
type  has  its  own  list  of  material  models.  For  the  2/D  and  3/D 
continuum  elements,  the  material  models  are  separated  into  two 
distinct  groups: 

1)  Linear  models  -  controlled  by  STSTL  (2/D  element)  or 
STST3L  (3/D  element). 

2)  Nonlinear  models  -  controlled  by  STSTN  (2/D  element) 
or  STST3N  (3/D  element). 

However,  for  truss  and  beam  elements,  coding  of  the  material  matrix 
are  interspaced  in  the  element  routines,  which  make  the  future 
program  changes  less  apparent. 

As  indicated  in  the  Appendix  C  of  ADINA  user  manual  [2], 
implementation  of  user-supplied  nonlinear  material  model  is  possible. 
Although  provisions  have  been  made  inside  the  program  to  accomodate 
new  material  model  for  2/D  and  3/D  elements,  such  implementation 
is  by  no  means  obvious.  The  difficulty  arises  from  the  dynamic 
allocation  of  working  storage  that  is  required  in  the  material 
subroutine.  A  clear  description  on  the  storage  allocation  was 
not  given  in  the  manual  and  therefore  the  user  must  figure  it 
out  from  the  existing  material  subroutines. 
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Since  the  constitutive  library  was  organized  as  a  lower  level 
module  of  the  element  library,  it  is  therefore  element-dependent. 

The  obvious  advantage  of  such  coding  structure  is  that  the  logical 
data  of  the  material  points  can  be  transferred  directly  to  the 
element  stiffness  calculations.  Moreover,  the  fast  core  storage 
for  material  data  processing  can  be  more  efficiently  utilized 
without  resorting  into  I/O  operations.  However,  this  advantage 
is  offset  by  the  following  shortcomings: 

1)  Codings  of  the  same  constitutive  model  for  different 
element  types  are  repetitive. 

2)  The  availability  of  constitutive  models  in  ADINA  is 
not  uniform  for  all  element  types. 

2.4  Output 

The  output  capability  of  ADINA  is  quite  simplistic,  it  only 
means  the  printing  of  nodal  displacements,  velocities  and  accelera¬ 
tions  (in  case  of  dynamic  analysis),  element  strains  and  stresses 
at  Gauss  points.  No  post-processing  capability  is  available, 
although  it  is  highly  desirable  for  the  large  volume  of  data 
generated  from  a  nonlinear  analysis. 

The  output  is  controlled  by  two  modules: 

WRITE  -  Printing  of  nodal  displacement,  velocities 
and  accelerations 

STRESS  -  Calculating  element  strains  and  stresses  at  Gauss 
points  and  printing  them  out. 

The  print-out  options  can  be  either  suppressed  or  made  at  designated 
time  blocks  according  to  the  input  instruction. 


A 
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3.  Use  of  Contro1  Variables 

For  a  given  structural  problem,  the  analysis  decision,  selection 
of  geometric  and  material  models  are  determined  in  ADINA  by  use  of 
conditional  statements  checking  the  values  of  control  variables.  If 
one  wishes  to  understand  the  macro  flow  of  the  program,  it  is  necessary 
to  know  the  precise  meaning  and  the  assigned  values  of  these  variables 
during  different  solution  stages. 

The  control  variables  used  in  ADINA  can  be  grouped  into  four  cate¬ 
gories: 

1)  General  Type  -  Variables  to  determine  solution  mode,  stress 
calculations  and  print-out,  time  history  count,  etc. 

2)  Analysis  Procedures  -  Variables  to  determine  static  or  dyna¬ 
mic,  linear  or  nonlinear  analysis  option. 

3)  Dynamic  Analysis  Detail  -  Variables  to  determine  the  type 
of  mass  system  and  integration  operator  to  be  used. 

4)  Stiffness  Reformation  &  Iterations  -  Variables  to  determine 
whether  the  stiffness  reformation  or  equilibrium  iteration 
is  desired. 

In  accordance  with  the  above  categories,  definitions  of  the  control 
variables,  subroutines  in  which  they  are  referenced,  are  listed  in 
Tables  11  -  14. 


Table  11.  Control  Variables  of  General  Type 


Name 

Definition 

Subroutine  Referenced 

KPRI 

Flag  fo>  stress  calculation  stage 
=  0  Yes,  calculate  element 
stresses^ 

=  1  No  stress  calculation 

AAMAIN 

TDFE2,  EL PAL 2 

KSTEP 

A  counter  to  determine  the  number 

of  load  (or  time)  step 

AAMAIN,  ASSEM,  BLKCNT 

KTR 

Flag  indicating  calculation  stage 
performed  in  equation  solver  "C0LS0L" 

*  1  Perform  matrix  triangulari- 

4 

zation  and  vector  solution 

*  2  Perform  vector  solution  only. 

AAMAIN  C0LS0L 

MODEX 

Flag  indicating  solution  mode 

*  0  Data  check  only 

*  1  Normal  execution 

=  2  A  restart  job 

AAMAIN,  AD  IN  I,  ASSEM, 

EQUIT,  LOADF.F,  NEWDAV, 

RSTART,  WRITE 

MODEL 

Flag  for  determining  material 
model  number 

=  1,  2,  ....  etc. 

ELCAL,  TODMFE2,  ST STL 3 

STSTN3,  INITWA3 ,  ELPAL 

NPAR1 

Flag  for  selecting  element  type 

=  1  Truss  element 

=  2  2/D  Continuum  element 

=  3  3/D  Continuum  element 

=  4  Beam  element 

ELCAL 

ELEMNT 

Notes:  1.  The  stress  calculation  option  is  superceded  by  an  input  variable 
IPS.  If  IPS  =  0,  no  stress  calculation  is  desired. 

2.  The  flag  is  also  used  in  other  equivalent  element  routines. 

3.  The  flag  is  also  used  in  other  equivalent  material  subroutines. 


4.  For  linear  analysis,  only  matrix  triangularization  is  performed. 
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Name 

IE  IG 

IND 


ISTAT 

KLIN 

INDNL 

I DEATH 


Table  12.  Control  Variables  for  Analysis 


Definiti in 


Flag  indicating  frequency  solution 

=  0  No  frequency  solution 
=  1  Yes 

Flag  to  control  various  stages  of 
calculations 

=  1  Form  linear  stiffness  matrix 
=  2  Form  mass,  damping  and  effec¬ 
tive  stiffness  matrices 
=  3  Perform  frequency  analysis 
=  4  Calculate  nonlinear  element 
stiffness  and  effective  load 
vector 

Flag  indicating  static  or  dynamic 
analysis 

=  0  Static  analysis 
=  1  Dynamic  analysis 

Flag  indicating  linear  or  nonlirear 
problem 

=  0  Linear  Problem 
=  1  Nonlinear  Problem 

Flag  indicating  linear  or  type  of 
nonlinear  analysis  (an  element 
level  flag) 

=  0  Linear  analysis 
=  1  Material  nonlinearity  with  small 
*  2  TL  formulation  large  ceformation 
=  3  UL  formulation  large  ceformation 

Flag  indicating  element  birth  and 
death  option  (an  element  group  flag) 

=  0  Elements  are  active 
a  1  Elements  have  birth  option 
=  2  Elements  have  death  option 


Procedures 

Subroutine  Referenced 
AAMAIN,  ADINI 

AAMAIN,  ADIDI ,  ASSEM, 

TODMFE1,  TDFE1,  QUADS1 
ELT2D32 

AAMAIN,  ADINI,  ASSEM, 
EQUIT,  INITAL,  LOADEF, 
NEWDAV,  RSTART ,  WRITE 

AAMAIN,  ADINI,  ASSEM, 
ELCAL,  LOADEF,  NEWDAV, 
RSTART 

ELCAL,  TDFE1 
QUADS1,  STSTN2 

deformation 

ELCAL,  TODMFE1,  TDFE1 


Notes:  1.  The  flag  is  also  used  in  other  equivalent  element  routines. 
2.  The  flag  is  also  used  in  other  material  subroutines. 
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Table  13.  Control  Variables  for  Dynamic  Analysis  Detail 


Name  Definition 


Subroutine  Referenced 


IMASS  Flag  indicating  type  of  mass  system  AAMAIN,  AD  INI ,  ASSEM, 

=  0  No  mass  effect  LOADEF,  TDFE1 

=  1  Lumped  mass  matrix  , 

=  2  Consistent  mass  matrix 


IDAMP  Flag  indicating  nodal  damping  effect  AD INI,  ASSEM 

=  0  No  damping  effect 
=  1  Yes 


IOPE 


Flag  indicating  time  integration 
operator 

=  0  Default  set  to  "2" 

=  1  Wilson  operator 

=  2  Newmark  operator 

=  3  Central  difference  operator 


AAMAIN,  ADINI,  ASSEM, 
INITAL,  LOADEF,  NEWDAV, 
RSTART. 


Note:  1.  Flag  is  also  used  in  other  equivalent  element  routines. 


■ 
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Table  14. 

Control  Variables  for  Stiffness  Reformation  and  Iterations 

Name 

Definition 

Subroutine  Referenced 

IREF 

Flag  indicating  reformation  of 
stiffness  matrixl 

=  0  Yes 
=  1  No 

AAMAIN,  ASSEM, 

BLKCNT,  TDFE2 

ICOUNT 

Flag  indicating  whether  calculation 
is  during  iteration  stage 

=  2  No 
=  3  Yes 

AAMAIN,  EQUIT,  TODMFE2, 
TDFE2,  QUADS2,  EL PAL 3 

IEQREF 

Flag  indicating  convergence  status 
in  iterations 

AAMAIN,  EQUIT 

=  0  Yes,  solution  has  converged 
*  1  No,  solution  does  not  converge. 

IEQUIT 

Flag  indicating  whether  iterations 
are  to  be  performed 

AAMAIN,  BLKCNT 

*  0  Yes 

=  1  No  iterations 

Notes:  1. 

Stiffness  reformation  may  >e  executed 
or  any  specified  number  of  time  steps, 
reassigned  in  subroutine  BLKCNT. 

every  time  (load)  step, 

The  value  of  IREF  is 

2. 

Flag  is  also  used  in  other  equivalent 

element  routines. 

3. 

Flag  is  also  used  in  other  material  subroutines. 

i 

t 


Commentary: 


1)  The  most  frequently  referenced  control  variables  are  "IND,  ISTAT, 
KLIN,  INDNL,  KSTEP,  MODEX,  MODEL".  These  a 'e  the  major  parameters 
to  determine  the  solution  phases  in  the  program. 

2)  Initial ization  of  the  control  variables  is  partly  assigned  in 
the  input  module  ADINI  and  partly  scattered  in  the  main  driver 
AAMAIN;  therefore  it  is  difficult  to  trace  the  exact  meaning  of 
these  variables.  For  example,  the  definition  of  the  analysis 
control  parameter  "IN0"  is  being  changed  intermittently  from  one 
solution  stage  to  the  other  and  the  reassignment  of  its  definition 
also  depends  on  the  status  of  control  variables  IE IG  and  MODEX. 

3)  As  seen  in  Table  12,  use  of  the  parameters  IE IG  vs.  ISTAT  and 
KLIN  vs.  INDNL  is  repetitive.  The  four  parameters  can  be  con¬ 
solidated  to  reduce  several  conditional  statements  in  the  code. 

4.  Data  Structure 

Typical  to  a  finite  element  software,  large  quantities  of  data 
are  being  processed  and  transferred  from  one  module  to  another  in 
ADINA.  This  is  done  by  use  of  both  fast  core  memory  and  secondary 
(low  speed)  storage.  The  fast  core  is  utilized  by  defining  a  dynami¬ 
cally  allocated  storage  array  in  the  blank  cormion  block  as  well  as 
fixed-dimension  arrays.  For  the  secondary  storage,  both  sequential 
files  and  random  (direct)  access  files  are  used.  No  clear-cut  data 
base  design  is  made  in  the  code,  nor  is  there  any  formal  data  manage¬ 
ment  system,  which  is  essential  for  large  scale  modern  softwares  [70,71]. 
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4.1  Fast  Core  Memory 

Two  types  of  In-core  array  allocations  were  made  in  ADINA : 

1)  Data  arrays  with  fixed  lengths  (or  dimensions) 

2)  A  variable  array  in  blank  common  block  which  is  randomly 
addressable. 

Only  a  small  portion  of  the  core  memory  is  blocked  into  fixed 
dimension  arrays.  Some  of  these  arrays  are  used  for  data  transfer 
through  labeled  common  blocks  between  lower  level  modules  whereas  other 
arrays  are  defined  locally  by  fixed  dimension  statements. 

Most  part  of  the  fast  core  is  effectively  utilized  by  defining 
a  master  array  in  the  blank  common  block  which  has  variable  length. 

This  array  is  then  partitioned  dynamically  into  subregions  for  the 
data  processing  and  data  storage  in  various  modules.  The  master  array 

is  primarily  used  as  the  work  space  and  therefore,  it  is  scratchable. 

Data  are  being  transferred  between  the  core  and  disk  or  tape  files 
independently  in  various  modules  without  going  through  a  data  manager. 
Moreover,  both  the  real  and  integer  numbers  are  stored  in  a  mixed  pattern. 
The  reason  for  using  dynamic  allocation  of  arrays  is  that  the  total  length 
and  the  manner  in  which  the  blank  common  block  is  partitioned  vary  with: 

1)  The  size  of  an  analysis  problem. 

2)  Type  of  analysis  option  requested  by  the  user. 

3)  Solution  stage  being  carried  out  in  the  program. 

Therefore,  it  is  difficult  to  follow  the  relative  position  and  status 
of  a  data  array  in  the  blank  common  block  unless  a  concise  storage 
picture  is  known.  For  this  reason,  the  dynamic  allocation  of  data 
arrays  in  the  blank  common  block  are  outlined  in  the  Appendix  B  ac¬ 
cording  to  the  following  solution  phases: 
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B1 .  Input  phase  -  nodal  data  and  initialization  of  array  variables. 
B2.  Input  phase  -  element  data. 

a.  Truss  element. 

b.  2/D  continuum  element. 

c.  3/D  continuum  element. 

d.  Beam  element. 

B3.  Addresses  of  diagonal  elements  of  structural  stiffness  matrix. 
B4.  Determination  of  blocks  for  structural  stiffness  matrix. 

B5.  Assembly  of  constant  structural  matrices,  i.e.  linear  stiff¬ 
ness,  mass  and  concentrated  damping  matrices. 

B6.  Construction  of  load  vectors. 

B7.  Data  for  a  restart  job. 

B8.  Calculation  of  effective  load  vector. 

B9.  Assembly  of  tangent  structural  stiffness  (for  out-of-core 
solution) . 

BIO.  Linear  equation  solver. 

Bll.  Equilibrium  iteration 

B12.  Static  analysis  or  dynamic  analysis  with  implicit-time 
integration. 

B13.  Dynamic  analysis  with  explicit  time  integration. 

The  block  diagrams  for  the  array  allocations  in  terms  of  the  posi¬ 


tion  of  the  pointer,  array  names  and  their  descriptions,  subroutines 
referenced  are  given  in  sections  B1  -  B13  of  Appendix  B. 
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Commentary: 

1)  Fast  core  memory  is  very  efficiently  utilized  by  dynamic  allocation 
of  arrays  in  the  blank  common  block. 

2)  No  central  data  management  is  available  for  data  transfer  between 
the  modules  or  between  the  fast  core  and  secondary  storage  devises. 
This  will  make  the  insertion  or  addition  of  new  arrays  into  the 
blank  common  block  much  more  difficult  without  over-writing  the 
existing  records. 

3)  The  efficiency  of  the  fast  core  usage  is  hampered  by  the  coding  of 
fixed  length  arrays  in  various  sub-modules  of  the  program. 

4.2  Low  Speed  Storage 

A  total  of  12  low  speed  storage  units  (not  including  the  input 
unit  #5  and  output  unit  #6)  are  used  in  ACINA.  Two  of  the  units  (#2 
and  #10)  are  random  (direct)  access  files  and  the  rest  of  the  units 
are  sequential  files.  Both  disk  and  tape  storage  units  are  utilized 
as  indicated  in  Table  15.  A  detailed  description  of  the  data  files 
for  storage  units  #1  -  #13  (not  including  #5  and  #6)  is  given  in 
12].  As  seen  from  Table  15,  most  of  the  data  files  are  of  the 

sequential  type,  which,  by  definition,  must  be  processed  sequentially. 
That  is,  data  records  must  be  searched  and  sorted  before  they  can  be 
processed.  Since  the  data  stored  on  these  units  are  fairly  straight¬ 
forward,  use  of  the  sequential  files  appears  to  be  quite  natural  for 
data  write  and  retrieval.  Nevertheless,  two  drawbacks  were  noticed: 
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Table  15. 

Disk/Tape  Utilization 

Unit 

Device 

Description 

1-13* 

Disk 

Data  files  (see  Table  16) 

56  \ 

Nodal  point  temperatures  prepro¬ 

59 

►  Tape 

cessed  input  data  porthole  files 

60 

for  saving  nodal/element  responses. 

63  J 

*  For 

a  restart  job. 

unit  #8  is  a  tape  file. 

J 
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1)  A  large  number  of  low  speed  storage  units  are  defined  and 
consequently  the  program  becomes  heavily  I/O  bound.  For 
example,  if  one  is  to  analyse  a  large  size  nonlinear  dynamic 
problem,  inordinate  amounts  of  data  transfer  between  memory 
and  out-of-core  files  are  necessary  which  will  in  turn  reduce 
the  computational  efficiency  of  the  code. 

2)  Since  a  module  or  sub-module  may  read  the  data  from  or  write 
the  data  to  the  sequential  files  independently,  any  inappro¬ 
priate  changes  in  the  data  structure  (dynamically  allocated) 
would  over-write  the  data  block  without  noticing  it. 

In  addition  to  the  use  of  sequential  files,  units  #2  and  #10  are 
the  random  access  files,  in  which  access  can  be  made  directly  to  any 
storage  location  by  computing  its  address  from  the  pre-determined  control 
number.  Unit  #2  stores  the  nonlinear  element  group  data,  such  as  ele¬ 
ment  connectivity  array,  material  constants,  strain  and  stress  com¬ 
ponents,  etc.  (see  Section  B2  );  certain  portions  of  these  data  are 
frequently  needed  during  the  various  stages  of  computations  whereas 
other  portions  (e.g.  element  stresses  and  strains)  have  to  be  updated. 
Moreover,  Unit  #10  stores  the  blocked  effective  stiffness  matrices 
from  which  considerable  effort  in  data  retrieval  is  necessary  due  to 
the  effect  of  block  coupling  during  equation  solving.  Therefore,  it 
is  also  a  natural  choice  to  make  both  units  #2  and  #10  as  the  random 
access  files  so  that  a  portion  of  the  records  on  the  file  can  be  con¬ 
veniently  processed  according  to  its  logical  sequence. 


V.  SUMMARY  AND  CONCLUSIONS 


The  finite  element  program  ADINA  was  reviewed  and  evaluated  with 
respect  to  its  nonlinear  analysis  capability  and  program  architecture. 

In  this  part  of  the  report,  discussion  was  focussed  on  the  theoretical 
basis,  numerical  approximations  adopted,  and  programing  details  from 
a  user's  point  of  view.  It  is  noted  that  this  evaluation  work  was 
made  on  the  1977  version  of  ADINA. 

Based  on  our  study,  some  of  the  unique  features  about  ADINA  can 
be  summarized: 

1)  "he  user's  manual  was  clearly  written  and  easy  to  follow.  Even 

ror  a  new  user,  the  input  cata  can  be  prepared  in  a  fairly  straight 
forward  manner  according  to  the  instructions  without  much  confusion 

2)  Theoretical  manual  is  well  documented,  in  which  the  definitions 
of  element  stiffness,  mass  and  material  matrices  used  in  the  code 
are  available  for  user's  who  wish  to  understand  the  code.  Most 
of  its  theoretical  development  can  be  found  from  the  developer's 
publications  referenced  in  [2]  and  [3]. 

3)  The  coding  practice  in  ADINA  corresponds  largely  to  the  theory 
outlined  in  its  manual. 

4)  Element  types  in  conjunction  with  material  models  were  organized 
on  a  modular  basis.  If  a  user  wishes  to  add  a  new  element  or 
material  model,  this  can  be  done  without  upsetting  other  parts 

of  the  program.  However,  addition  of  new  elements  must  be  limited 
to  those  having  6  or  less  than  6  degrees  of  freedom  at  each  node. 
Otherwise,  modification  of  the  code  becomes  less  apparent. 
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5)  Calculations  of  element  stiffness  matrices,  especially  for  simple 
material  models  such  as  linearly  elastic,  isotropic  materials, 
were  coded  by  explicit  matrix  multiplications.  In  this  way,  the 
amount  of  computational  effort  is  kept  minimum. 

6)  The  fast  core  memory  is  efficiently  utilized  by  dynamic  allocations 
of  the  blank  common  block,  which  has  a  variable  length  depending 

on  the  size  of  the  problem  and  type  of  analysis  requested  by  the 
user. 

7)  The  program  is  relatively  machine  independent  (primarily  for  main¬ 
frames).  Double  precision  arithemetic  option  for  the  IBM  or  UNIVAC 

•  computers,  and  necessary  coding  changes  for  various  machines  are 
included  in  the  program  by  commented  statements. 

8)  Organization  of  the  program  and  data  transfer  between  the  modules 
were  designed  primarily  according  to  the  computational  flow  of 
the  analysis  procedures.  Therefore,  ADINA  is  a  high  effictncy 
code-  for  executions. 

9)  With  the  out-of-core  linear  equation  solver,  ADINA  can  be  used 
for  the  analysis  of  large  size  application  problems.  For  this 
purpose,  the  code  runs  more  efficiently  on  a  virtual  core  environ¬ 
ment. 

10)  Since  ADINA  is  a  relatively  "small"  general  purpose  finite  element 
program  (e.g.  compared  to  NASTRAN)  and  many  commented  statements 
are  included  in  the  program,  it  is  well  suited  as  a  research  tool 
for  studies  in  structural  mechanics. 
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Nevertheless,  AOINA  has  its  limitations.  The  limitations  result 

largely  from  the  narrow  scope  for  which  the  program  was  developed. 

A  summary  of  these  limitations  are  given  as  below: 

1)  The  greatest  drawback  of  the  program  is  the  lack  of  a  comprehensive 
pre-  and  post-  processing  package.  Especially  the  preprocessing 
capability  in  terms  of  grid  generation,  grid  plot  and  thorough 
data  check  is  essential  for  the  nonlinear  analysis  of  application 
problems. 

2)  The  1977  version  of  ADINA  has  basically  four  different  element 
types:  3/D  truss,  2/D  isoparametric  continuum,  3/D  isoparametric 
continuum,  and  3/D  beam.  Therefore,  geometric  representation 

of  the  code  for  complex  structures  is  somewhat  limited. 

3)  The  manner  in  handling  large  deformation  problems  is  not  complete. 

In  fact,  the  current  version  can  only  be  used  for  large  displacements 
but  small  strains. 

4)  The  incremental  solution  technique  adopted  in  ADINA  is  a  modified 
Newton-Raphson  method,  for  which  the  tangent  modulus  together 
with  equilibrium  iterations  is  used.  This  method  fails  to  handle 
the  situation  where  the  load-deformation  response  of  a  structure 
exhibits  softening  behavior  and  then  suddenly  becomes  stiffened. 

5)  When  ADINA  is  applied  to  nonlinear  dynamic  problems,  there  is 
always  a  question  in  determining  the  size  of  the  time  increment 
for  obtaining  convergent  and  accurate  solutions.  For  application 
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purpose,  it  is  necessary  to  include  a  heuristic  self-adaptive 
scheme  so  that  the  program  can  determine  the  convergent  solution 
increments  internally  without  user's  guesswork. 

5)  Coding  organization  of  the  procedure  library  does  not  conform 
to  the  modular  concept.  Consequently,  it  is  rather  difficult 
to  alter  this  part  of  the  coding. 

7)  Data  transfer  between  the  modules  or  between  the  fast  core  and 
secondary  storage  units  is  made  without  a  structured  data  manage¬ 
ment.  Therefore,  it  is  difficult  to  follow  the  status  of  the 
data  files  or  to  make  any  changes  on  these  data. 

Due  to  the  aforementioned  limitaions,  ADINA  is  considered  primarily 
a  research  tool  for  nonlinear  finite  dement  structural  analysis.  Although 
additional  elements,  material  models  or  other  features  can  be  added 
to  the  program,  it  is  doubtful  that  the  code  could  be  expanded  signifi¬ 
cant!/  beyond  its  current  level. 
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Appendix  A 


FLOW  DIAGRAMS  FOR  VARIOUS  SOLUTION  PHASES 


(Analysis  phase)  1  (Input  phase) 


A1 


A1  Flow  Diagram  of  the  Main  Program  "AAMAIN" 
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2.  Flow  Diagram  for  Reading  Input  Data 
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A3.  Processing  of  Element  Data 
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A7.  Equilibrium  Iteration 
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A9.  TRUSS  ELEMENT 


A10.  2-D  Solid  Element 
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Appendix  B 


DYNAMIC  ALLOCATIONS  OF  ARRAYS  IN  THE 
BLANK  COMMON  BLOCK 


B1.  Input  Phase  -  Nodal  Data  and  Initialization  of  Arrays 


Calling  Subroutine:  AD INI 


Address  Dimension  INPUT  NODMAS  INITAL 


Description 


Nodal  point  Coordinates: 


NUMNP*NDOF  ID 


Boundary  condition  array 


NUMNP* ITWO  X 


Nodal  point  X-coordinate 


NUMNP* I TWO  Y 


Nodal  point  Y-coordinate 


NUMNP* ITWO 


Nodal  point  Z-coordinate 


Initial  Conditions: 


NUMNP*NDOF 


Boundary  condition  array 


NEQ*ITW0 


DISP  Displacement  input 


NEQ*ITW0 


Velocity  input 


NEQ*ITW0 


Acceleration  input 


NEQ*ITW0 


DISPM2  l.  Nodal  mass  vector  or  Vodal 
damping  vector 

2.  For  explicit  time  integration, 
the  displacement  vector  at 


N1  =  NBCST  +2*NEGNl  +2 

NBCST  =  Variable  to  shift  start  of  BLANK  COMMON  beyond  last  address  of  program  storage. 
NEGNL  p  Number  of  nonlinear  element  group. 


Address 


NUMNP*ITWO 

NUMNP* I TWO 
‘NUMNP*ITWO 

NEQ+1 


Input 

Plase  - 

Elen ent 

Data 

'outines 

:  TRUSS 

,  TODMFE 

,  THREDM, 

BEAM 

TRUSS 

TODMFE 

THREDM 

BEAM 

Description 

ID 

ID 

ID 

ID 

Boundary  condition  array 

X 

- 

X 

X 

Nodal  point  X-coordinate 

Y 

Y 

Y 

Y 

Nodal  point  Y-ct ordinate 

Z 

Z 

Z 

Z 

Nodal  point  Z-ccordinate 

HT 

HT 

HT 

HT 

Column  heights 

TEMPV1 

TEMPV1 

TEMPV1 

TEMPV1 

Nodal  point  temperatures 

EE 

EE 

EE 

EE 

Element  group  information 

MAXEST 


Note:  If  temperature  is  not  involved.  No  =  N5  +  NEQ. 


B2-A.  Element  Data  for  Truss 


Address 


NFIRST 


NT  01 
=  N102 


Nl  03 

N105 

N106 


N107 

N108 

N109 

N110 

Nlll 

Nil  2 

N113 

Nil  4 


N115 
N1 16 


Calling  Subroutine:  TRUSS 


Dimension 


20 


NUMMAT*ITWO 


NUMMAT*ITWO 


NDM*NUME 


NDM*NUME*ITWO 


NUME 


NUME* I TWO 


NUME 


MM*NUME*ITWO 


MM*NDM*  NUME* I TWO 


NUME* I DWA* I TWO 


NCON*NUMMAT*ITWO 


NDW*NUME*MXNODS 


NUME 


RUSS 


NPAR 


DEN 


AREA 


LM 


XYZ 


MATP 


EPSIN 


IPS 


PROP 


■{ 


Description 


Element  group  control  parameters 


Mass  density 


Cross-sectional  area 


Element  connectivity  array  (vector) 


Element  nodal  coordinates 


Material  property  set  number 


Initial  strain 


Stress  printing  flag 


ETIMV  Element  birth/death  time 


EDISB  Element  disps.  at  birth  time 


WA=IWA  Working  array 


Material  constants 


NODGL  Global  node  numbers 


IELTD  ELEMENT  number  of  nodes 


N6  INDj*  4 


First  Address:  .1FIRST  =  U10  JN|V  4  ,  MM  ID£ATH  =  x 


[0  I  DEATH  =  0 


Last  Address:  NLAST  =  N116 
Total  Length  of  Record:  MIDEST  =  NLAST-NFIRST  +  1 


B4 


B2-B.  Element  for  2-D  Element 


Address 


Calling  Subroutine:  TOOMFt 
Dimension  TDFE 


Description 


NFIRST 

N101 

N102 

N103 

N104 

N105 

N106 

N107 

N108 

N109 

N 1 1 0 
Nlll 

N112 

N113 

N114 


NPAR 


LM 
Y  L 

1ELT 

I  PST 
BETA 

THICK 

MATP 

DEN 

PROP 

WA 

NOD5 

ETIMV 

EDISB 
I  TABLE 


Element  croup  control  parameters 

Element  connectivity  array 
Element  nodal  coordinates 

Element  number  of  nodes 

Stress  printing  flag 
Material  angle 

Element  thickness 

Material  property  set  number 
Mass  density 

Material  constants 

Working  array 

Midside  nodes  location  array 

Element  expiry  time  array 

Element  birthtime  nodal  coordinates 
Stress  output  location  tables 


NLAST 


B2-C.  Element  Data  for  3-D  Element 


Calling  Subroutine: 

Address  Dimension  THDFE 

NPAR 
LM 

XYZ 
IELTD 
IELTX 
IPST 

MATP 

N0D9 
1  REUSE 
DEN 
PROP 

WA 
ETIMV 

EDISB 
ITABLE 

DCA 

MAXES V 


THREDM 

Description 

Element  group  control  parameters 
Element  connectivity  array 

Element  nodal  coordinates 

Element  number  of  nodes  for  snape 
function 

Number  of  nodes  used  to  describe 
elemerc  geometry 

Stress  printing  flag 

Material  property  set  number 

Midside  nodes  location  array 
Index  for  repeating  element 

Mass  density 

Material  constants 

Working  array 

Element  expiry  time  array 

Element  birthtime  nodal 
coordi nates 

Stress  output  location  tables 

Direction  cosine  array 

Material  axis  orientation 
storage  vector 


NFIRST 

N101 

N102 

N103 

N104 

N105 

N107 

N108 

N109 

N110 

Nlll 

N112 

N113 

N114 

N115 

N116 
N 1 17 

NLAST 


B2-C.  Element  Data  for  3-D  Element  (Continued) 


Cal  lino  Subroutine  THREDM 


Address  Dimension  THDFE 


S 

XM 

B 

RE 

EDIS 


Description 


Element  stiffness  matrix 
Nodel  mass 

Compacted  strain-disp.  matrix 
Element  nodal  force  vector 
Element  displacement  vector 


B7 


B2-D.  Element  Data  for  BEAM 


Calling  Subroutine:  BEAM 


Address  Oimensioi 


BMEL 


Description 


NPaR 

E 

G 

DEN 

XI 

Y I 

ZI 

AREA 

XYZ 

LM 

IPS 

MATP 

PRUP 

ICS 


Element  group  control  parameters 
Young's  modulus 
Shear  modulus 
Mass  density 

Second  moment  of  area  @  '  R  axi  ; 

Second  moment  of  area  @  'S'  axis 

Second  moment  of  area  @  ' T*  axis 

Normal  +  Shear  section  area 

Element  coordinates  array 
Element  connectivity  array 
Stress  output  flag 

Material  property  set  number 
Material  constants 
Section  identification  flag 


N114 


NFAC*NUMMAT 


1SHEAR 


Flag  for  transverse  shear  effects 


B2-D.  Element  Data  for  BEAM  (Continued) 


Cal  lino  Subroutine:  BEAM 


Address  Dimension  BMEL 


Description 


WA 

ITABLE 

SR 

ETIMV 

EDISP 

IMOMNT 

IELRET 

PDISP 

GAMA 


Working  array 
Stress  output 

Gauss  elimination  coefficient 
Element  expiry  array 
Element  disps.  at  birth  time 
Number  of  end  release  tables 
Element  number  of  end  release  table 
Element  displacement  vector 
Creates  element  transformation  matrix 


% 


! 


B9 


B3.  Addresses  of  Diagonal  Elements  of  Stricture  1  Stiffness  Matrix 

Calling  Subroutine:  ADDRES 

Address  Dimension  Address  Description 


Locations  of  diagonals  of 
structural  stiffness  matrix 


MHT 


Col u nn  height 


Note:  If  temperature  is  involved,  N6  =  N5  +  NEQ  +  1. 


B4.  Determination  of  Blocks  for  Stiffness  Matrix 


Calling  Subroutine:  SBLOCK 


Address  Dimension 


SBLOCK 

Description 

MAXA 

Locations  of  diagonal 

of 

structural  stiffness  i 

natrix 

NCOLBV 

Number  of  columns  per 

block 

ICOPL 

First  coupling  block 

BIO 


B5.  Assembly  of  Constant  Structural  Matrices 


Calling  Subroutine:  ASSEM 


Address  Dimension 


ASSEM 


Description 


MAX  A 

Locations  of  diagonals  of 
structural  stiffness  matrix 

NCOLBV 

Number  of  columns  per  block 

ICOPL 

First  coupling  blocks 

IGRBLC 

Nonlinear  element  group  to 
stiffness  Mock  coupling 

For  the  case  of  frequency  analysis 

AA~^ 

^•Storing  blocks  of  matrices 

CC^ 

DD 

Concentrated  masses  and  lumped 
mass  matrix 

BB 

Concentrated  dampers 

EE 

Element  group  storage 

(Storage  for  element  calculations) 


Notes: 


1.  N1  =  NBCST+2*NEGNL*2 

2.  For  one  block  case,  N1D  =  NIC 

3.  For  implicit  time  integration  aid  one  block  case  N4  =  N3 
For  explicit  time  integration  N4  =  N1 

4.  For  static  analysis  N6  =  N4 

For  static  analysis  with  lumped  gravity  load  N6  =  N5. 


l. 


Bll 


B6.  Construction  of  Load  Vectors 


Calling  Subroutines:  LOAD,  CLOADS,  GRAVL,  TODPRL,  THDPRL 


Address  Dimension 


Arrays 


description 


MAXA 

ID 

RG 

RGST 

R 


Locations  of  diagonal:,  of 
structural  stiffness  matrix 

Boundary  condition  array 

Interpolated  values  of  time 
tunctiois 

Initial  values  of  time  functions 
Load  vector 


Notes : 

M4  plus 

(2*NPTM  +  2*NL0AD  +  NEQ)*ITW0 

in  calculation 

of  concentrated 

+  (NEQ  +  4*NL0AD) 

load  vector 

M4  plus 

(3*NUMNP  +  NDFR2*NPR2  +  4*NPR2)*ITW0 

in  calculation 

of  2/D  pressure 

+  NPR2*(NDFR2  +  6) 

load  vector 

M4  plus  (3*NUMNP  +  NDFR3  +  5*NPR3)*ITWO 
+  NPR3*(NDFR3  +  N0DE3  +  2) 


in  calculation  of  3/D  pressure 
load  vector 


B7.  Transfer  of  Data  for  a  Restart  Job 


Address 


Calling  Subroutine:  RSTART 
Dimension  RSTART 


Description 


N2  - - 

NEQ*ITW0 

N7  - - — 

NEQ*ITW0 

N8  - - — 

"*  NEQ*ITW0 

N10 - 

"*  MAXEST 

DISP 

VhL 

ACC 

EE 


Initial  displacement 
Initial  velocity 
Initial  acceleration 
Element  information 


B13 


Address 


B8.  Calculation  of  Effective  Load  Vector 
Calling  Subroutine:  LOADEF 

Dimension  LOADEF  Description 


MAXA 

NCOLBV 


DISP 

R 

AA 


WV 


VEL 

ACC 

XM 


Locations  of  diagonals  of 
structural  stiffness  matrix 

Number  of  columns  per  block 


Initial  displacement 
Load  vector 

Linear  part  of  structural 
stiffness 


Working  vector 


Nodal  velocity 
Nodal  acceleraton 
Nodal  mass 


Notes:  1.  For  one  block  case,  N1D=N1C,  N4B  =  N4A 

2.  For  explicit  time  integration  N2  =  N1  +  NEQMTWO,  N5  =  N4 


B14 


B9.  Assembly  of  Tangent  Structural  Stiffness 
Calling  Subroutine:  ASSEM 


Address  Dimension  ASSEM 


Description 


MAXA 

NCOLBV 


Locations  of  diagonals  of 
structural  stiffness  matrix 

Number  of  columns  per  block 


IGRBLC  Nonlinear  element  group  to 

stiffness  block  coupling 


DD 

BB 


Concentrated  masses  and  lumped 
mass  matrix 

Concentrated  dampers 


Storing  blocks  of  matrices 


TEMPV2  Nodal  point  temperatures 


EE  Element  group  storage 


i 


BIO.  Linear  Equation  Solver 


Calling  Subroutine  -  COLSOL 


Dimension  COLSOL 


Description 


MAXA 

NCOLBV 

ICOPL 


V 

A 

B 

D 


Locations  of  diagonals  of 
structural  stiffness  matrix 

Number  of  columns  per  block 

First  coupling  block 


Load  or  displacement  vector 
Structural  stiffness  matrix 

Working  space  for  blocked 
structural  stiffness 

Reduced  diagonal  location 


I 


\  1 


B16 


Bll.  Equilibrium  Iteration 
Callinq  Subroutine:  EQUIT 


Address  Dimension  EQUIT 


Description 


N1 

N1A- 

NIB 

NIC" 

N1D' 

N2 

N3- 

N4 

N4A 

N4B ' 

N5  ■ 


NEQ  +  1 


NBLOCK 


NBLOCK 


NBLOCK*NEGLNL 


( IEIG+l)*NBLOCK+l 


NEQ*ITWO 


NEQ*ITWO 


ISTOH*ITWO 


ISTOH*ITWO 


NEQ*ITW0 


NEQ*ITWO 


MAXA 

NCOLBV 

ICOPL 


DISP 

DISPI 

AA 

CC 

DK 

RE 


Addresses  of  diaconal  elements  in 
effective  stiffness  matrix 

Number  of  columns  per  block 


First  coupling  block 


Displacement  at  previous  time  step 


Displacement  increment  at  previous 
time  step 

Effective  stiffness  matrix  and 
working  storage 

Working  storage  in  out-of-core 
solution 

Elements  of  D  in  LDLT  factorization 
of  effective  stiffness  matrix 

Out  of  balance  loads 


i 


i 


t 


I 

B16 

(continued) 

Bll.  Equilibrium  Iteration 
Calling  Subroutine:  EQUIT 

Address  Dimension  EQUIT  Description 


Working  Array 

Velocity  at  previous  time  step 
Acceleration  at  previous  time  step 
Nodal  mass 

Element  group  storage 


B17 


B12.  Storage  Allocation  During  Static  Analysis  OR 
Implicit-  Time  Integration 


Address  Dimension 


N1 

N1A 

NIB 

NIC 

N2 

N3 

N4 

N4A 

N4B 

N5 

N6 

N6A 

N6B 

N7 

N8 

N9 

N10 


NEQ+1 


NBLOCK 


NBLOCK 
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Arrays 

MAXA 

NCOLBV 

ICOPL 

IGRBLC 

DISP 

R 

A 

B 

Dk 

RE 

WV 

TEMP VI 
TEMP V 2 
VEL 
ACC 
XM 
EE 


Description 

Locations  of  diagonals  of  structural 
stiffness  matrix 

Number  of  columns  per  block 

First  coupling  block 

Nonlinear  element  group  to  stiffness 
block  coupling 

Displacement  vector 

Load  vector  and  displacement  increments 


Storage  for  solution  of  equations 


Our-of  balance  loads,  and  correction  to 
displacement  increment 

Working  vector 


Working  vectors  to  share  nodal  point 
temperatures 


Velocity  vector  for  dynamic  analysis  only 


Acceleration  vector  for  dynamic  analysis 
only 

Lumped  mass  matrix  for  dynamic  analysis 
only 

Element  group  information 


Total 


B18 


B 1 3 .  Storage  All( nation  During  Explicit  Time  Integration 


Address  Dimension 


Arrays 

DISPM 

DISP 

R 

W  V 


TEMPV1* 


TEMPV2- 


VEL 

ACC 

EE 


Description 

Displacement  vector  at  time  t  -  At 

Displacement  vector  at  time  t 

Load  vector  at  time  /  or  displacement 
vector  at  time  t  +  At 

Working  vector 

Working  vectors  to  store  nodal  point 
temperatures 

Velocity  vector  at  time  t 
Acceleration  vector  at  time  t 
Element  group  information 
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